# Full Delay Logistic Population Model with Sustainable Harvesting

Authors: Martin Anokye; Elvis Kobina Donkoh; Michael Fosu Ofori; Agnes Adom-Konadu
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

The paper studies the dynamics of a full delay-logistic population model incorporated with a proportionate harvesting function. The study discusses the stability of the model in comparison with the well-known Hutchinson logistic growth equation with harvesting function using the rate of harvesting as a bifurcation parameter to determine sustainable harvesting rate even at a bigger time delay ofτ=3.00. In all cases, the Hutchinson equation with harvesting was forced to converge to equilibrium using an additional and a different time delay parameter, a deficiency previous researchers have failed to address when the Hutchinson model is used for this purpose. The population fluctuations are catered for with this model making the estimated maximum sustainable growth and harvest reflect realities as this model drastically reduces time-delay associated oscillations compared to the well-known Hutchinson delayed logistic models. The numerical simulations were be done using the MatLab Software.


---

## Body

## 1. Introduction

Logistic growth model has been considered in different contexts among researchers [1–4]. The simplest form of a logistic model to study population dynamics is given by
(1)ẋt=rxt1−xtc,forxt0=x0.Thext is considered as the size of the population at time t, where the r parameter represents the intrinsic growth rate, and the positive constant c is the systems’ carrying capacity i.e., the maximum population the system can sustain [5]. Recent population studies have witnessed a surged interest in dynamics of population retarded models involving delay logistic equations. Many attempts have also been made to find reasonable mathematical equations with time lags to describe complex physical systems [6–9].For example, Hutchinson in [2] used equation (2) to address the aforementioned issue incorporating time delay in logistic population model (1). In the preceding model, the per capita rate of change of the system depends on the state of the system in the past.
(2)ẋt=rxt1−xt−τc,fort>0,(3)xt=φt,fort∈−τ,0,where r,τ,c are positive numbers, with xt=φt, where φt∈−τ,0,ℝ as the initial history function. Another time delay parameter was introduced in equation (2) to obtain the following logistic growth equation with the same initial condition as stated in equation (2),
(4)ẋt=rxt−τ1−xt−τc,where r the growth rate, c the population or environmental capacity, and τ time delay are all positive constants. This model was studied in [3, 10, 11]. Usually, in most populations, outside influences affect the proportion of lives within the system. Such influences may be caused by factors including harvesting, drugs for treatment, and poaching just to mention a few. Assuming population modeled by equation (1) is affected at a rate proportional to its size; then, we get a different model defined by
(5)ẋt=rxt1−xtc−hxt,where h≥0 is the effect of the influence on population at time t. This gives the fraction of the population that is harvested [12]. Similarly, Berezansky et al. in [13] incorporated h in equation (2) and obtained the following equation
(6)ẋt=rxt1−xt−τc−hxt−τ,t∈0,τ.In this model, the rate of outside influence is proportional to the population size at timet−τ. This model supports the arguments that any interacting species in the system are dependent on amassed resources and hunting or harvesting effects in the past.Similarly, Piotrowska and Bodnar in [3] and Cooke et al. [14] used the model below by introducing time delay in rxt to postulate that the intrinsic growth rate depends on past time t−τ, where τ, the developmental time of the population, is in the system with model given by
(7)ẋt=rxt−τ1−xt−τc−htxt,t∈0,τ,where h is constant function at a step time t. The difference between equations (6) and (7) is that the rate of proportionality of outsides interference on the population is not dependent on the time in the past t−τ.From the two equations, we can infer that no author has yet considered time delay in all parts of the model as given by(8)ẋt=rxt−τ1−xt−τc−hxt−τ,t∈0,τ,where the population growth rate, the per capita rate of change, and the effect of the outside influence are all determined on time in the past.In this study, we will consider equation (8) with time delays in all parts. We will also explore the possible effects of changing the model parameters as means to compare equations (8) and (6) and then discuss the defects in the Hutchinson equation (incorporated with outside interference) which previous authors have so far failed to address when the model is used with harvesting function. Berezansky et al. in [13] studied the existence, positivity, and the sufficient conditions for the extinction of the solution of the equation (6), but this study will beyond this by estimating the sustainable growth and harvest rate. This is the first of its kind that a research study has touched on the deficiencies of the well-known Hutchinson equation, especially when using it as a harvesting incorporated function. The objective of the study is to help the research community with the selection of the right model for the right purpose. We will identify the dynamical properties of equation (8) so that managers can adopt it as an interventional plan to manage the population. Thus, we will use harvest rate h as the bifurcation parameter for the system. Since h is regarded as the rate of outside interference; then, ecological managers can use it in making decisions on harvesting efforts. The aforementioned bifurcation strategy is adopted from the argument by Geritz and Kisdi [15] that parameters from population models are helpful in determining the state of the models.The paper is structured as follows: Section2 discusses the existence, stability, and bifurcation of equation (8). The section would also consider conditions for positivity of the solution in equation (8). In Section 3, detailed numerical solutions and stability analyses would be done comparing the models (8) and (6) by applying the same parameter values. The numerical simulations will be done using the MatLab Software. Sections 4 outlines the findings and conclusion derived from the analyses of the models.


## 2. Materials and Methods

### 2.1. Existence, Stability, and Bifurcation of Equation (8)

As we normalize the system, then, we can observe that ifh<r, the two possible equilibrium of equation (8) will be xt=0 and xt=1−h/r. The second part xt=1−h/r is stable and so for any initial condition x0>0, all solutions of equation (8) will converge towards it as t⟶∞. We let ch denote the positive equilibrium of the equation (8) to emphasize the fact that h is the bifurcation parameter of the system as stated in the objective of this study. On the other hand, if h≥r, then, all solutions will be attracted to the former equilibrium point x=0 as the only solution of equation (8). This means the population should be driven towards extinction.We now transform equation (8) into the following form by the change of variables using Xt=r/r−hxt(9)Ẋt=−r−hXt−τ−1+Xt−τc,fort>0.From equation (9), the linearized form is expressed as follows;
(10)Ẋt=−r−hXt−τ.Denotingxt=κeγt as solution to equation (10) can result in the following transcendental equation whose zero solutions are negative real parts and is supposed to occur at xt=cr−h/r (i.e., the only positive equilibrium supposed to be ch=1−h/r);
(11)γ+ρe−γτ=0,where ρ=r−h. Setting γ=μ+iδ, for μ∈ℝ, δ∈0,∞, then, we have the following pairwise equation.
(12)μ=−ρe−μτcosδτ,δ=ρe−μτsinδτ.From equation (11), we obtain a negative real eigenvalue γ=−ρ<0 if we let τ=0. However, if τ increases such that the real part γ=μ>0, then, for some τ0, there results in the real part γτ0=μτ0=0. This provides a pair of imaginary roots ±iδ0,δ0=δτ0 (i.e., from equation (11)), which confirms the assertion in [2]. We can infer that cosδ0τ=0 and consequently obtain δ0τk=π/2+2kπ,k=0,1,2,⋯Sinceδ0=ρ=r−h, then by extension, we can conclude that the only positive equilibrium applicable to equation (9) is ch=1−h/r, and it is stable if r−hτ<π/2, and unstable when r−hτ>π/2. This implies that for rτ>π/2 with which nonharvested model (4) is unstable, then becomes stable if the rate of harvesting reaches the critical value h∗=r−π/2τ.


#### 2.1.1. Stability Switches and Hopf Bifurcation of (8)

Also assumingρ>0 (i.e., r−h>0); then, stability switches are determined from the complex eigenvalues γ=iδ where (δ>0 is a solution of equation (11)). Since sinδτ is supposed to be positive (refer to equation (12)); then, with ρ=δ, we have the stability switches given by
(13)δρ=π2+2kπk=0,1,2,⋯.Consideringγ as function of the bifurcation parameter ρ, we can find the direction of the stability switches by differentiating equation (11) with respect ρ as expressed implicitly by
(14)γ̇+e−γτ−ρτe−γτγ̇=0.This implies that(15)γ̇=e−γττρe−γτ−1=γρτγ+ρ=iδiδρτ+ρ.Therefore(16)γ̇=iδρ−iδρτρ+iδρτρ−iδρτ=iδρ−iδρτδτρ2+ρ2,where the real part δ2τρ/δτρ2+ρ2>0 shows that the eigenvalue has changed from being negative to a positive real number [2].At the emergence of Hopf bifurcation, we letμ=0, so that γ=iδ, then, from equation (11), we obtain
(17)iδ=−ρe−iδτ=−ρcosδτ−isinδτ.Equating the real and imaginary parts of this equation results in the following(18)0=ρcosδτ,(19)δ=ρsinδτ.By squaring both sides of the two adding them together, we have(20)δ=∓ρ.The set of solutions ofτ associated with these equations is given by
(21)τn=π/2+2nπδ,forδ>0,3π/2+2nπδ,forδ<0,where n=0.1.2,⋯. If we let δ0=ρ, then, γnτ=μnτ+iδnτ is the root of the equation (11) satisfying μn=0 and δnτn=δ0.We now have the following results.Lemma 1.
μ̇nτn>0.Proof.
In differentiating both sides of equation (11) with respect to τ, we have
(22)dγdτ=γρe−γτ−ρe−γτ=−ρe−γτ−dγdτ−γ.Given thatγ=−ρe−γτ,(23)dγdτ=−γ21+γτ,(24)dγdττ=τn=δ021+iδ0τ=δ02−iδ03τ1+iδ02τ.Therefore, this implies that(25)μ̇nτn=δ021+δ02τ>0.This results confirms that, in lettingγ=μ+iδ, if μ varies from the left to right, then, for μ<0, the system is in a stable state, μ>0 a bifurcation state, and μ=0 a limiting case [2, 16].


#### 2.1.2. Existence and Uniqueness Solution of (8)

From equation (9), we let a=r−h, c=1; then, we obtain the following initial value problem
(26)Ẋt=−aXt−τ−1+Xt−τ,fort>0,(27)Xt=φt,fort∈−τ,0,where a,τ are positive numbers, with Xt=φt, as the initial history function. We consider the interval 0,τ to generate a nonnegative solution to equation (26) from its equivalent expression given by
(28)Xt=φ0+−a∫0tXs−τ−1+Xs−τds.Sinceφ0≥ and a>0, the solution exists and it is unique and nonnegative in the neighbourhood considered. Again in the interval t∈n−1τ,nτ, if we let Xn:n−1τ,nτ⟶ℝ+ be the solution of equation (26); then, for t∈nτ,n+1τ, it implies that
(29)Xt=Xnnτ+−a∫nτtXns−τ−1+Xs−τds.Thus, it is observed that for every nonnegative initial functionφt≥0, the solution of equation (26) is defined for t≥0 in t∈nτ,n+1τ.


#### 2.1.3. Existence and Positivity Solution of (8)

Following the Theorem 1.2 in [10], the solution of (26) can have negative values for a positive initial condition. Therefore, due to biological constraints, we study the conditions which will guarantee nonnegative solutions of equation (26) for every positive initial function φt. If we let t=t¯τ and then drop the bar, then after transformation of the equation (26), we obtain the following initial value problem
(30)Ẋt=−aτXt−1−1+Xt−1,fort>0,(31)Xt=φt,fort∈−1,0.As a consequence of the Theorem 1.2 in [9] as stated earlier, we adopt the following.Let(32)0≤φt≤1,fort∈−1,0,due to biological reasons and(33)φt=12,fort∈−1,0,1,fort=0,then on the intervaln−1,n, we let xn be the solution of equation (30). We can observe that
(34)X1t=φ0+∫0t−aτφs−1+φsds.Hence,(35)X1t=1+∫0taτ4ds=1+aτt4,(36)X2t=X11+∫0t−1−aτX1s−1+X1sds=1+aτ4+∫0t−1−aτ1+aτs4aτs4ds.Therefore, we obtain(37)X2=X22=1+aτ4−aτ28−aτ348=P1aτ.If the condition (32) is satisfied then following inequality also holds:
(38)∀t>0,Xt≤1+aτ4.IfXt>1, then, there exists t0<t such that Xt0=1, and
(39)Xt=Xt0+∫t0−1t−1−aτXs−1+Xsds,so that we can have(40)Xt≤1+∫t0−1t0aτ4ds=1+aτ4,or else we have
(41)Xt≥1+∫t0−1t0−aτ1+aτ4aτ4ds=P2sds=1−aτ24−aτ316=P2aτ.Then, as a consequence of Theorem 2.1 in [10], if aτ<r2, it implies the polynomial P2aτ>0 and P1aτ<0 for aτ>r1, where r1 and r2 are the greatest roots of P1 and P2, respectively. That is r1=3.0578 and r2=1.6786.


## 2.1. Existence, Stability, and Bifurcation of Equation (8)

As we normalize the system, then, we can observe that ifh<r, the two possible equilibrium of equation (8) will be xt=0 and xt=1−h/r. The second part xt=1−h/r is stable and so for any initial condition x0>0, all solutions of equation (8) will converge towards it as t⟶∞. We let ch denote the positive equilibrium of the equation (8) to emphasize the fact that h is the bifurcation parameter of the system as stated in the objective of this study. On the other hand, if h≥r, then, all solutions will be attracted to the former equilibrium point x=0 as the only solution of equation (8). This means the population should be driven towards extinction.We now transform equation (8) into the following form by the change of variables using Xt=r/r−hxt(9)Ẋt=−r−hXt−τ−1+Xt−τc,fort>0.From equation (9), the linearized form is expressed as follows;
(10)Ẋt=−r−hXt−τ.Denotingxt=κeγt as solution to equation (10) can result in the following transcendental equation whose zero solutions are negative real parts and is supposed to occur at xt=cr−h/r (i.e., the only positive equilibrium supposed to be ch=1−h/r);
(11)γ+ρe−γτ=0,where ρ=r−h. Setting γ=μ+iδ, for μ∈ℝ, δ∈0,∞, then, we have the following pairwise equation.
(12)μ=−ρe−μτcosδτ,δ=ρe−μτsinδτ.From equation (11), we obtain a negative real eigenvalue γ=−ρ<0 if we let τ=0. However, if τ increases such that the real part γ=μ>0, then, for some τ0, there results in the real part γτ0=μτ0=0. This provides a pair of imaginary roots ±iδ0,δ0=δτ0 (i.e., from equation (11)), which confirms the assertion in [2]. We can infer that cosδ0τ=0 and consequently obtain δ0τk=π/2+2kπ,k=0,1,2,⋯Sinceδ0=ρ=r−h, then by extension, we can conclude that the only positive equilibrium applicable to equation (9) is ch=1−h/r, and it is stable if r−hτ<π/2, and unstable when r−hτ>π/2. This implies that for rτ>π/2 with which nonharvested model (4) is unstable, then becomes stable if the rate of harvesting reaches the critical value h∗=r−π/2τ.


### 2.1.1. Stability Switches and Hopf Bifurcation of (8)

Also assumingρ>0 (i.e., r−h>0); then, stability switches are determined from the complex eigenvalues γ=iδ where (δ>0 is a solution of equation (11)). Since sinδτ is supposed to be positive (refer to equation (12)); then, with ρ=δ, we have the stability switches given by
(13)δρ=π2+2kπk=0,1,2,⋯.Consideringγ as function of the bifurcation parameter ρ, we can find the direction of the stability switches by differentiating equation (11) with respect ρ as expressed implicitly by
(14)γ̇+e−γτ−ρτe−γτγ̇=0.This implies that(15)γ̇=e−γττρe−γτ−1=γρτγ+ρ=iδiδρτ+ρ.Therefore(16)γ̇=iδρ−iδρτρ+iδρτρ−iδρτ=iδρ−iδρτδτρ2+ρ2,where the real part δ2τρ/δτρ2+ρ2>0 shows that the eigenvalue has changed from being negative to a positive real number [2].At the emergence of Hopf bifurcation, we letμ=0, so that γ=iδ, then, from equation (11), we obtain
(17)iδ=−ρe−iδτ=−ρcosδτ−isinδτ.Equating the real and imaginary parts of this equation results in the following(18)0=ρcosδτ,(19)δ=ρsinδτ.By squaring both sides of the two adding them together, we have(20)δ=∓ρ.The set of solutions ofτ associated with these equations is given by
(21)τn=π/2+2nπδ,forδ>0,3π/2+2nπδ,forδ<0,where n=0.1.2,⋯. If we let δ0=ρ, then, γnτ=μnτ+iδnτ is the root of the equation (11) satisfying μn=0 and δnτn=δ0.We now have the following results.Lemma 1.
μ̇nτn>0.Proof.
In differentiating both sides of equation (11) with respect to τ, we have
(22)dγdτ=γρe−γτ−ρe−γτ=−ρe−γτ−dγdτ−γ.Given thatγ=−ρe−γτ,(23)dγdτ=−γ21+γτ,(24)dγdττ=τn=δ021+iδ0τ=δ02−iδ03τ1+iδ02τ.Therefore, this implies that(25)μ̇nτn=δ021+δ02τ>0.This results confirms that, in lettingγ=μ+iδ, if μ varies from the left to right, then, for μ<0, the system is in a stable state, μ>0 a bifurcation state, and μ=0 a limiting case [2, 16].


### 2.1.2. Existence and Uniqueness Solution of (8)

From equation (9), we let a=r−h, c=1; then, we obtain the following initial value problem
(26)Ẋt=−aXt−τ−1+Xt−τ,fort>0,(27)Xt=φt,fort∈−τ,0,where a,τ are positive numbers, with Xt=φt, as the initial history function. We consider the interval 0,τ to generate a nonnegative solution to equation (26) from its equivalent expression given by
(28)Xt=φ0+−a∫0tXs−τ−1+Xs−τds.Sinceφ0≥ and a>0, the solution exists and it is unique and nonnegative in the neighbourhood considered. Again in the interval t∈n−1τ,nτ, if we let Xn:n−1τ,nτ⟶ℝ+ be the solution of equation (26); then, for t∈nτ,n+1τ, it implies that
(29)Xt=Xnnτ+−a∫nτtXns−τ−1+Xs−τds.Thus, it is observed that for every nonnegative initial functionφt≥0, the solution of equation (26) is defined for t≥0 in t∈nτ,n+1τ.


### 2.1.3. Existence and Positivity Solution of (8)

Following the Theorem 1.2 in [10], the solution of (26) can have negative values for a positive initial condition. Therefore, due to biological constraints, we study the conditions which will guarantee nonnegative solutions of equation (26) for every positive initial function φt. If we let t=t¯τ and then drop the bar, then after transformation of the equation (26), we obtain the following initial value problem
(30)Ẋt=−aτXt−1−1+Xt−1,fort>0,(31)Xt=φt,fort∈−1,0.As a consequence of the Theorem 1.2 in [9] as stated earlier, we adopt the following.Let(32)0≤φt≤1,fort∈−1,0,due to biological reasons and(33)φt=12,fort∈−1,0,1,fort=0,then on the intervaln−1,n, we let xn be the solution of equation (30). We can observe that
(34)X1t=φ0+∫0t−aτφs−1+φsds.Hence,(35)X1t=1+∫0taτ4ds=1+aτt4,(36)X2t=X11+∫0t−1−aτX1s−1+X1sds=1+aτ4+∫0t−1−aτ1+aτs4aτs4ds.Therefore, we obtain(37)X2=X22=1+aτ4−aτ28−aτ348=P1aτ.If the condition (32) is satisfied then following inequality also holds:
(38)∀t>0,Xt≤1+aτ4.IfXt>1, then, there exists t0<t such that Xt0=1, and
(39)Xt=Xt0+∫t0−1t−1−aτXs−1+Xsds,so that we can have(40)Xt≤1+∫t0−1t0aτ4ds=1+aτ4,or else we have
(41)Xt≥1+∫t0−1t0−aτ1+aτ4aτ4ds=P2sds=1−aτ24−aτ316=P2aτ.Then, as a consequence of Theorem 2.1 in [10], if aτ<r2, it implies the polynomial P2aτ>0 and P1aτ<0 for aτ>r1, where r1 and r2 are the greatest roots of P1 and P2, respectively. That is r1=3.0578 and r2=1.6786.


## 2.1.1. Stability Switches and Hopf Bifurcation of (8)

Also assumingρ>0 (i.e., r−h>0); then, stability switches are determined from the complex eigenvalues γ=iδ where (δ>0 is a solution of equation (11)). Since sinδτ is supposed to be positive (refer to equation (12)); then, with ρ=δ, we have the stability switches given by
(13)δρ=π2+2kπk=0,1,2,⋯.Consideringγ as function of the bifurcation parameter ρ, we can find the direction of the stability switches by differentiating equation (11) with respect ρ as expressed implicitly by
(14)γ̇+e−γτ−ρτe−γτγ̇=0.This implies that(15)γ̇=e−γττρe−γτ−1=γρτγ+ρ=iδiδρτ+ρ.Therefore(16)γ̇=iδρ−iδρτρ+iδρτρ−iδρτ=iδρ−iδρτδτρ2+ρ2,where the real part δ2τρ/δτρ2+ρ2>0 shows that the eigenvalue has changed from being negative to a positive real number [2].At the emergence of Hopf bifurcation, we letμ=0, so that γ=iδ, then, from equation (11), we obtain
(17)iδ=−ρe−iδτ=−ρcosδτ−isinδτ.Equating the real and imaginary parts of this equation results in the following(18)0=ρcosδτ,(19)δ=ρsinδτ.By squaring both sides of the two adding them together, we have(20)δ=∓ρ.The set of solutions ofτ associated with these equations is given by
(21)τn=π/2+2nπδ,forδ>0,3π/2+2nπδ,forδ<0,where n=0.1.2,⋯. If we let δ0=ρ, then, γnτ=μnτ+iδnτ is the root of the equation (11) satisfying μn=0 and δnτn=δ0.We now have the following results.Lemma 1.
μ̇nτn>0.Proof.
In differentiating both sides of equation (11) with respect to τ, we have
(22)dγdτ=γρe−γτ−ρe−γτ=−ρe−γτ−dγdτ−γ.Given thatγ=−ρe−γτ,(23)dγdτ=−γ21+γτ,(24)dγdττ=τn=δ021+iδ0τ=δ02−iδ03τ1+iδ02τ.Therefore, this implies that(25)μ̇nτn=δ021+δ02τ>0.This results confirms that, in lettingγ=μ+iδ, if μ varies from the left to right, then, for μ<0, the system is in a stable state, μ>0 a bifurcation state, and μ=0 a limiting case [2, 16].


## 2.1.2. Existence and Uniqueness Solution of (8)

From equation (9), we let a=r−h, c=1; then, we obtain the following initial value problem
(26)Ẋt=−aXt−τ−1+Xt−τ,fort>0,(27)Xt=φt,fort∈−τ,0,where a,τ are positive numbers, with Xt=φt, as the initial history function. We consider the interval 0,τ to generate a nonnegative solution to equation (26) from its equivalent expression given by
(28)Xt=φ0+−a∫0tXs−τ−1+Xs−τds.Sinceφ0≥ and a>0, the solution exists and it is unique and nonnegative in the neighbourhood considered. Again in the interval t∈n−1τ,nτ, if we let Xn:n−1τ,nτ⟶ℝ+ be the solution of equation (26); then, for t∈nτ,n+1τ, it implies that
(29)Xt=Xnnτ+−a∫nτtXns−τ−1+Xs−τds.Thus, it is observed that for every nonnegative initial functionφt≥0, the solution of equation (26) is defined for t≥0 in t∈nτ,n+1τ.


## 2.1.3. Existence and Positivity Solution of (8)

Following the Theorem 1.2 in [10], the solution of (26) can have negative values for a positive initial condition. Therefore, due to biological constraints, we study the conditions which will guarantee nonnegative solutions of equation (26) for every positive initial function φt. If we let t=t¯τ and then drop the bar, then after transformation of the equation (26), we obtain the following initial value problem
(30)Ẋt=−aτXt−1−1+Xt−1,fort>0,(31)Xt=φt,fort∈−1,0.As a consequence of the Theorem 1.2 in [9] as stated earlier, we adopt the following.Let(32)0≤φt≤1,fort∈−1,0,due to biological reasons and(33)φt=12,fort∈−1,0,1,fort=0,then on the intervaln−1,n, we let xn be the solution of equation (30). We can observe that
(34)X1t=φ0+∫0t−aτφs−1+φsds.Hence,(35)X1t=1+∫0taτ4ds=1+aτt4,(36)X2t=X11+∫0t−1−aτX1s−1+X1sds=1+aτ4+∫0t−1−aτ1+aτs4aτs4ds.Therefore, we obtain(37)X2=X22=1+aτ4−aτ28−aτ348=P1aτ.If the condition (32) is satisfied then following inequality also holds:
(38)∀t>0,Xt≤1+aτ4.IfXt>1, then, there exists t0<t such that Xt0=1, and
(39)Xt=Xt0+∫t0−1t−1−aτXs−1+Xsds,so that we can have(40)Xt≤1+∫t0−1t0aτ4ds=1+aτ4,or else we have
(41)Xt≥1+∫t0−1t0−aτ1+aτ4aτ4ds=P2sds=1−aτ24−aτ316=P2aτ.Then, as a consequence of Theorem 2.1 in [10], if aτ<r2, it implies the polynomial P2aτ>0 and P1aτ<0 for aτ>r1, where r1 and r2 are the greatest roots of P1 and P2, respectively. That is r1=3.0578 and r2=1.6786.


## 3. Results and Discussion

### 3.1. Numerical Solution: Stability Analysis of Model (2) and (4)

We now do numerical simulations with the models (2) and (4) (from which we derived the equations (6) and (8), respectively) to examine their dynamics regarding the effect of the delay parameters incorporated in them. Berezansky et al. in [13] studied the existence, positivity, and the sufficient conditions for the extinction of the solution of the equation (6), but this study will beyond this by estimating the sustainable growth and harvest rate. We will compare the two equations so that we will be able to recommend the one to select when using the model as a harvest incorporated model.Forτr=1.56>π/2, derived from τ=11.00,r=0.15, as shown Figure 1, the system exhibited periodic solution at x∗=c=1. The results of the simulations from the two plots for the given parameter values complement Theorem 2.1 in [2].Figure 1
Stability analysis of the model (2).


#### 3.1.1. Numerical Solution: Stability Analysis of Model (4)

Since equation (4) is further fused with another delay parameter, we want to study if it adheres to the conditions in Theorem 2.1 in [2], as done with equation (2).In Figure2, we used the same values for Figure 1, as displayed in the legends. This time, the space occupied by the various oscillations is reduced as shown by the population axis. This confirms that when the time delay is applied in all parts of the model, the oscillations associated with delays are rather minimized.Figure 2
Stability analysis of the model (4).


### 3.2. Numerical Solution: Harvesting Parameter for Bifurcation

We study the stability dynamics of equations (8) and (6) using the same parameters together with the bifurcation parameter h. In this section, we will use the MatLab software for all our numerical simulations. To normalize the system, we let x0=0.5, and c=1, then, with r=0.1, from equation (5), we have
(42)ẋt=xt0.11−xt−τ−hxt−τ.At equilibrium,xt=xt−τ=xe, and there exists two equilibrium solutions:
(43)xe=0,andxe=1−10h.Usingxe=1−10h, we can determine the positive equilibrium of the system since xe is the fish population, and it should not be negative. Assuming the harvest rate h>0.1 and xe=1−10h is negative; then, the harvest will drive the fish population to extinction. On the other hand, if h=0.1 (same as the growth rate r); then, the results will affirm the other equilibrium solution xe=0, and all the fish population will disappear. The equilibrium solutions are also applicable to equation (8) as shown in equation (9) using the parameters r−hτ and π/2 and also as we let xe=ch=1−10h. The figures below illustrate this assertion.In Figure3, for r=h=0.1, the result of equation (8) conforms to that of the xe=0, and therefore, all the total fish population will disappear with time but not same with equation (6). Equation (6) moves away from the zero equilibrium and can only be forced to it when we increase the time delay associated with rate of harvest or introduce a different time delay on the harvest rate.Figure 3
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.1.In Figure4, a new time delay was introduced on the harvesting effort so that equation (6) has two different discrete time delays and now we have both equation (6) and equation (8) converging towards the zero equilibrium point xe=0, which denotes that the total fish population with time will disappear.Figure 4
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.1.Similarly, if the harvest rate is set ath=0.025, then, xe=ch=1−100.025=0.75, which is stable. It is also observed that with 0≤h<0.1, the equilibrium fish population will always be positive and can be predicted and they will not die. The simulation of the harvest rate in the above-given interval is demonstrated in Figures 5 and 6.Figure 5
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.025.Figure 6
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.025.From Figure5 above, it is observed that for a proportional decrease in fish population through harvesting by h=0.025, there is a substantial growth of 0.75 as shown in the the plot. However, the model (6) moves away higher above the point, and it has to be redirected.In Figure6, it can be concluded that as the harvesting rate decreases, positive equilibrium is established, and there is an upsurge in the average and variations in population size. This time, both models converged to the equilibrium point xe=0.75 because equation (6) was forced to do so.


### 3.3. Optimal Sustainable Harvesting Value

Now, given the fish population sizexe=ch, the total harvest can be defined by
(44)hch=h1−10h.We will determine the maximum sustainable harvest at equilibrium at the level that optimizesh1−10h. Let Wh denote the total function of harvest rate h, then, from equation (44), we have
(45)Wh=h1−10h.We compute the local maximum of the equation (45) as follows; Ẇh=1−20h, then, we have h=0.05 as the critical point. We also take the second derivative to find the type of point. W¨h=−20<0, which denotes the local maximum point. This implies that at h=0.05, as a fraction of the fish population offers the maximum sustainable harvest value. Therefore, with fish growth rate r=0.1 and optimal harvest rate h=0.05, we obtain ch=1−100.05=0.5, which is exactly half of the maximum carrying fish population capacity. This is proven in Figure 7 that follows.Figure 7
Stability analysis of equation (8) with r=0.1, τ=3.00, and h=0.05.In Figure7, it is indicated that if we set the harvest rate at h=0.05, and the initial condition at x0=0.5, then, the system will provide its maximum sustainable growth and harvest. This happens at exactly half of the maximum carrying fish population capacity. The two models settled at the same maximum level because equation (6) has been forced to do so through the other different time delay incorporated in the model.In Figure8, it is also observed that if the harvest rate is set at h=0.05, with the initial condition at x0=1.0, then, the system will converge to exactly half the maximum carrying fish population capacity which denotes the level within which we can harvest. The two models have settled at the same level because equation (6) has been adjusted with another time delay parameter.Figure 8
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.05.


## 3.1. Numerical Solution: Stability Analysis of Model (2) and (4)

We now do numerical simulations with the models (2) and (4) (from which we derived the equations (6) and (8), respectively) to examine their dynamics regarding the effect of the delay parameters incorporated in them. Berezansky et al. in [13] studied the existence, positivity, and the sufficient conditions for the extinction of the solution of the equation (6), but this study will beyond this by estimating the sustainable growth and harvest rate. We will compare the two equations so that we will be able to recommend the one to select when using the model as a harvest incorporated model.Forτr=1.56>π/2, derived from τ=11.00,r=0.15, as shown Figure 1, the system exhibited periodic solution at x∗=c=1. The results of the simulations from the two plots for the given parameter values complement Theorem 2.1 in [2].Figure 1
Stability analysis of the model (2).


### 3.1.1. Numerical Solution: Stability Analysis of Model (4)

Since equation (4) is further fused with another delay parameter, we want to study if it adheres to the conditions in Theorem 2.1 in [2], as done with equation (2).In Figure2, we used the same values for Figure 1, as displayed in the legends. This time, the space occupied by the various oscillations is reduced as shown by the population axis. This confirms that when the time delay is applied in all parts of the model, the oscillations associated with delays are rather minimized.Figure 2
Stability analysis of the model (4).


## 3.1.1. Numerical Solution: Stability Analysis of Model (4)

Since equation (4) is further fused with another delay parameter, we want to study if it adheres to the conditions in Theorem 2.1 in [2], as done with equation (2).In Figure2, we used the same values for Figure 1, as displayed in the legends. This time, the space occupied by the various oscillations is reduced as shown by the population axis. This confirms that when the time delay is applied in all parts of the model, the oscillations associated with delays are rather minimized.Figure 2
Stability analysis of the model (4).


## 3.2. Numerical Solution: Harvesting Parameter for Bifurcation

We study the stability dynamics of equations (8) and (6) using the same parameters together with the bifurcation parameter h. In this section, we will use the MatLab software for all our numerical simulations. To normalize the system, we let x0=0.5, and c=1, then, with r=0.1, from equation (5), we have
(42)ẋt=xt0.11−xt−τ−hxt−τ.At equilibrium,xt=xt−τ=xe, and there exists two equilibrium solutions:
(43)xe=0,andxe=1−10h.Usingxe=1−10h, we can determine the positive equilibrium of the system since xe is the fish population, and it should not be negative. Assuming the harvest rate h>0.1 and xe=1−10h is negative; then, the harvest will drive the fish population to extinction. On the other hand, if h=0.1 (same as the growth rate r); then, the results will affirm the other equilibrium solution xe=0, and all the fish population will disappear. The equilibrium solutions are also applicable to equation (8) as shown in equation (9) using the parameters r−hτ and π/2 and also as we let xe=ch=1−10h. The figures below illustrate this assertion.In Figure3, for r=h=0.1, the result of equation (8) conforms to that of the xe=0, and therefore, all the total fish population will disappear with time but not same with equation (6). Equation (6) moves away from the zero equilibrium and can only be forced to it when we increase the time delay associated with rate of harvest or introduce a different time delay on the harvest rate.Figure 3
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.1.In Figure4, a new time delay was introduced on the harvesting effort so that equation (6) has two different discrete time delays and now we have both equation (6) and equation (8) converging towards the zero equilibrium point xe=0, which denotes that the total fish population with time will disappear.Figure 4
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.1.Similarly, if the harvest rate is set ath=0.025, then, xe=ch=1−100.025=0.75, which is stable. It is also observed that with 0≤h<0.1, the equilibrium fish population will always be positive and can be predicted and they will not die. The simulation of the harvest rate in the above-given interval is demonstrated in Figures 5 and 6.Figure 5
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.025.Figure 6
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.025.From Figure5 above, it is observed that for a proportional decrease in fish population through harvesting by h=0.025, there is a substantial growth of 0.75 as shown in the the plot. However, the model (6) moves away higher above the point, and it has to be redirected.In Figure6, it can be concluded that as the harvesting rate decreases, positive equilibrium is established, and there is an upsurge in the average and variations in population size. This time, both models converged to the equilibrium point xe=0.75 because equation (6) was forced to do so.


## 3.3. Optimal Sustainable Harvesting Value

Now, given the fish population sizexe=ch, the total harvest can be defined by
(44)hch=h1−10h.We will determine the maximum sustainable harvest at equilibrium at the level that optimizesh1−10h. Let Wh denote the total function of harvest rate h, then, from equation (44), we have
(45)Wh=h1−10h.We compute the local maximum of the equation (45) as follows; Ẇh=1−20h, then, we have h=0.05 as the critical point. We also take the second derivative to find the type of point. W¨h=−20<0, which denotes the local maximum point. This implies that at h=0.05, as a fraction of the fish population offers the maximum sustainable harvest value. Therefore, with fish growth rate r=0.1 and optimal harvest rate h=0.05, we obtain ch=1−100.05=0.5, which is exactly half of the maximum carrying fish population capacity. This is proven in Figure 7 that follows.Figure 7
Stability analysis of equation (8) with r=0.1, τ=3.00, and h=0.05.In Figure7, it is indicated that if we set the harvest rate at h=0.05, and the initial condition at x0=0.5, then, the system will provide its maximum sustainable growth and harvest. This happens at exactly half of the maximum carrying fish population capacity. The two models settled at the same maximum level because equation (6) has been forced to do so through the other different time delay incorporated in the model.In Figure8, it is also observed that if the harvest rate is set at h=0.05, with the initial condition at x0=1.0, then, the system will converge to exactly half the maximum carrying fish population capacity which denotes the level within which we can harvest. The two models have settled at the same level because equation (6) has been adjusted with another time delay parameter.Figure 8
Stability analysis of equations (6) and (8) with r=0.1, τ=3.00, and h=0.05.


## 4. Conclusion

The paper studies the dynamics of a full delay logistic population model with sustainable harvesting that is proportional to the fish population. In this study, we used the harvest rate as a bifurcation parameter to determine the maximum permissible harvest per fishing period. We compared the stability dynamics of our model with the well-known Hutchinson equation with the harvesting function. The model (4) (without harvesting) is found to be less oscillatory compared to the Hutchinson equation (2) as shown in Figure 2. In all cases, the Hutchinson equation with harvesting was found to be deficient as it had to be adjusted (incorporate a different delay parameter on the harvesting rate) before it could converge to an equilibrium point set to the system. Since the logistic growth model without time delay cannot account for fluctuations in future populations, it is believed that the estimated maximum sustainable harvest value obtained with equation (8) reflects the realities of fish population dynamics. The paper recommends that further studies are done considering the full delayed logistic population model, especially when using it as a harvesting incorporated function.



---

*Source: 6777827-2022-08-26.xml*