Spatio-temporal solutions of a diffusive directed dynamics model with harvesting.

The study considers a directed dynamics reaction-diffusion competition model to study the density of evolution for a single species population with harvesting effect in a heterogeneous environment, where all functions are spatially distributed in time series. The dispersal dynamics describe the growth of the species, which is distributed according to the resource function with no-flux boundary conditions. The analysis investigates the existence, positivity, persistence, and stability of solutions for both time-periodic and spatial functions. The carrying capacity and the distribution function are either arbitrary or proportional. It is observed that if harvesting exceeds the growth rate, then eventually, the population drops down to extinction. Several numerical examples are considered to support the theoretical results. Supplementary Information: The online version contains supplementary material available at 10.1007/s12190-022-01742-x.

Introduction

A classical reaction-diffusion equation (RDE) arises in many areas of sciences and engineering, such as chemistry [1], pattern formation [2, 3], ecological invasions [4], spread of epidemics [5], biophysics [6, 7], and many other areas [8-10]. In general, an RDE represents a mathematical model where several components interact. In ecology, a classical RDE model plays an important role in studying different types of species dynamics [10].

The biological diversity of the earth is affected by habitat loss and fragmentation. Understanding the consequences of these effects on the persistence of native populations and communities is one of the challenging problems in conservation biology. Interactions of two species may be altered due to habitat loss and fragmentation, which illustrates the critical links between spatial ecology and conservation [11]. For instance, during an avian influenza outbreak, different culling strategies could be incorporated as a control measure to minimize the long-term effects on genetic diversity [12]. Environmental factors do not only affect to the biological diversity but they influence human diseases too. Eco-epidemiological models are used to investigate the transmission dynamics of an infectious disease engaged with a prey-predator connection in the impact of sickness in ecological models [13-15]. A model consists of a four-dimensional system of differential equations is considered to study the global dynamics for a honeybee colony infested by virus-carrying Varroa mites in [16].

In theoretical ecology, the determination of distribution (density) of species (or organisms) and the structure of communities is one of the main goals where species (or organisms) interact with each other in a given environment [17]. In the recent decade, the role of spatial effects in preserving biodiversity has received a great deal of attention in the literature on conservation [18-21].

There is a substantial difference between classical and directed diffusion. For instance, in classical diffusion some movements are unpredictable; whereas a directed diffusion focuses on selecting habitat influences. Populations differ in the diffusion strategies they employ as well as in their environmental intensities [22-24]. Further, the competition between populations also has influences on the diffusion [25-27].

Both autonomous and non-autonomous population models subject to either impulsive or continuous harvesting were studied in [28, 29] and found that the impulsive strategy can be as good as the continuous one but cannot outperform. The effect of periodic habitat fluctuations on a nonlinear insect population model has been studied in [30] and illustrate the possibility of a significant increase in population numbers using both theorically and experimentally. The harvesting problem in the seasonal environment was considered in [31], presented a general harvesting model with seasonal carrying capacity and developed a theoretical foundation. A diffusive predator-prey models are analyzed and studied in [32, 33], complex dynamics of these models are visualized in the presence of the predator harvesting. A density-dependent prey-predator model having the effect of harvesting is investigated in [34]. Moreover, prey-predator population model with fear and Allee effects is considered in [35]. This paper investigates the global existence of the solution of a diffusive directed dynamic (DDD) model, which is a variant of the classical RDE model. In the proposed DDD model, we consider the resource distribution is uniform per capita according to the distribution function with the homogeneous Neumann boundary condition which is interpreted as there is no flux of population on the boundary; or the number of populations that come in and go out at each point on the boundary is equal at each moment, so the total number of populations is in balance. In addition, we incorporate the harvesting effect in the DDD model and study its dynamics.

In the present DDD model, we consider a dispersed function-driven diffusion when species migrate to its favorable environment, not just the carrying capacity. Note that the desired environment of any species can be equivalent to the carrying capacity, but the converse is not true [36]. There can be any partition of carrying capacity where the resource for any population is high. Thus, species move to that location from where it exists. The stability of solutions of the equation with the logistic type of growth was first studied in [37] and then extended in [38] with general growth function. The influence of harvesting for a two-species competition was studied in [39] with a spatially heterogeneous environment. In this paper, we consider the time and space involving disperse diffusion.

One important observation is how diffusion terms are added to an RDE or a DDD model. Commonly, in a regular diffusion, it takes the form , where d is the specific diffusion rate and v(t, x) is the density of the species at time t and position x. In the case of directed diffusion, carrying capacity-driven diffusion is an important term, and has the form , where K(t, x) is the carrying capacity of the environment. In this paper we incorporated the new directed diffusion term, , known as the disperse diffusion, where P(t, x) is some positive and bounded resource distribution function. Both K, and P have the same dimension as v.

Many organisms, which follow this type of diffusion strategy, compete for common land, food habitat, etc., slow each other’s growth rate and diffuse according to the properties of the common space they need. For instance, like deers, they tend to move toward green areas for the food chain. If a forest area is in-habitat, they move to other spaces or become extinct. This study explores a population distribution that has a tremendous application in different areas of ecology and economy; instantaneous examples are in river and ocean ecology to observe the seasonal behaviors of different species, including wild animals and winter birds. We investigate global well-posedness of the mathematical model, determine conditions on harvesting rate for which non-trivial equilibrium states exist and examine their global stability for both spatial and time dependent periodic functions with the same period. For time periodic growth, resource, capacity and harvesting functions, we prove the existence of time-periodic states.

The main objectives of the study are: (1) We will study with an influencing diffusion model where the resource distribution is uniform per capita instead of random to observe that the coherent strategy is favorable for the species in competition. (2) It is expected that any species’ movement (diffusion) according to the distribution function changes the persistence and extinction threshold. (3) Including the new diffusion strategy in the model leads the species to be depending on the resource distribution function, which governs the environmental carrying capacity. (4) In addition, theory does not give any idea about the profile of nonzero equilibrium state , except when resource function P is proportional to the capacity function K. We have explored this numerically since it is important from ecological considerations.

The rest of the paper is organized as follows: In Sect. 2, the DDD model is developed. In Sect. 3, the existence and uniqueness of solution of a DDD model is presented along with a convergence analysis. In Sect. 4, we present the proof of the existence of time-periodic solution in a bounded region. The dynamics of harvesting and its applications are being described in Sect. 5. Some examples are presented in Sect. 6 for applications. Finally, Sect. 7 concludes the summary and discussion of the results.

Mathematical model

To investigate the spatial effects in ecology, the random RDE models are well established and used widely (see [1]–[6]). These models support the following three important types of ecological phenomena: (1) A minimal patch size is needed to sustain a population and its growth, (2) Biological invasions generate the propagation of wavefronts, and (3) The spatial patterns are formed in the distributions of populations in homogeneous environments. By observing the classical diffusion model, in this paper, we consider the following DDD system to model the density of the species v(t, x):where is an open domain, , n is the outward normal vector to the boundary. P(t, x) is a resource distribution function, K(t, x) is the carrying capacity of the environment, r(t, x) is the specific growth rate, is the harvesting function, is the initial density of the species, for any , and . In this model, the distribution function ensures the resource is uniformly distributed per capita, also the population is assumed to be distributed according to the available per resource. The modified diffusion term leads to the homogeneous Neumann boundary condition, and was first introduced in [18] motivated by the choice of optimal harvesting strategies. This alternative model of time-dependent spatial reaction-diffusion equation to study the population distribution in a heterogeneous environment, where diffusive transport of population is proportional to the gradient of population density per unit capacity instead of just the population density. In this model, populations are distributed based on available resources. It is an influencing diffusion model where the resource distribution is uniform per capita to observe that the coherent strategy is favorable for the species in competition. Other nonlinear diffusion types, based on the experimental data, were considered in population dynamics [40-42].

Remark 1

The functions P, K, and r are all depend on both t, and x, unless otherwise state.

Throughout this paper, we assume in , for some and functions P, , K, and are uniformly bounded from below and above in .

The equilibrium point of the above system is , where is the obvious as a trivial steady-state. Our primary interest is non-trivial equilibrium. In absence of harvesting (), all possible equilibrium points of the above system can be obtained as either for or where k is a constant. If , then is the only non-trivial equilibrium point. If the resource function and capacity function are not proportional then is the non-trivial steady-state. For non-zero harvesting (), is the unique steady-state neither K(x) nor P(x). It is also noted that when diffusion coefficient is small, equilibrium state , whereas when diffusion effects are large, is proportional to P(x) as presented numerically.

We define , and and have the following assumptions: It is easy to verify that the logistic growth function with positive bounded r follows the above (h1)–(h3) assumptions. In the following sections, for establishing several results, we assume the conditions (h1)-(h3) hold.

h is uniformly Hölder continuous in , where J is either a suitable sub-interval of or a sector between lower and upper solutions. Moreover, functions h and are as functions of t, Hölder continuous in in the x variable and continuous in as functions of v.

g(t, x, v) is strictly monotonically decreasing in v, when , and when , for .

.

Existence, uniqueness, and convergence analysis

For convenience, everywhere in this section, we assume The domain is nonempty open bounded with . We first state and proof the existence and uniqueness theorem for the time dependent distribution function, carrying capacity, and specific growth rate.

P and K are Hölder continuous in x and continuously differentiable periodic function, for any ,

r is continuous in , periodic and for any .

Theorem 1

Let the initial function , , and in a nonempty open bounded sub-domain . Then there exists a unique attractor v of the problem (1)-(3) and the solution is positive.

Proof

For simplicity and to translate the diffusion term in classical style, we first make the following substitutionand since the function P is positive and bounded from above in , the new variable w is well defined. We obtain that function w satisfiesThen the equivalent system of (1)-(3) iswhere we defineFor some particular functions in the domain, either K is proportional to P or they are non-proportional. It is obvious that the operator is uniformly elliptic with the coefficient that is Hölder continuous. Therefore, according to Lemma A2, to show that the existence of the unique solution of (4)-(6), we only need to construct an ordered pair of upper and lower solutions. To construct an upper solution, denoteand consider a constant C such thatand which can be found according to (h2). Note that since is bounded in and K is bounded from below,Next, chooseand let . Thengivesfor some particular bounded functions P, and K if condition (7) holds. Thus, the first inequality of (57) in Definition A2 holds. Further,Therefore, is an upper solution of (4)-(6) by Definition A2. The function is obviously a lower solution. Check that the right-hand side of the following functionis continuously differentiable with respect to w and we can denote the maximal derivative of l in w for each (t, x)Then, the Lipschitz condition (59) holds and by Lemma A2, there exists a unique solution of the problem (4)-(6) satisfying . Finally, using the inverse substitution we obtain a unique solution of (1)-(3).

To ensure that the solution is positive for any non-negative initial function , we assumewhereThe inverse substitution, into (4)-(6), we haveTo study (11)-(13), note that s is non-negative since . Then, assume that s has a zero at some point . If then according to the maximum principle, we obtain in , for details, see Lemma A3. According to the assumption of the theorem on some open bounded sub-domain and since w is a solution of (4)-(6). Therefore in . Moreover, if for some then by Lemma A2, we have at which contradicts the no flux boundary conditions in (11)-(13). This concludes the proof.

Now we consider the system (1)-(3) without harvesting and with spatially distributed parameters P(x), K(x), and r(x):The result of Theorem 1 will apparently valid for (14)-(16). The next result shows that when the carrying capacity and the intrinsic growth rate are only distributed spatially, and then is an equilibrium solution of (14)-(16). Also, the solution v(t, x) of (1)-(3) with converges to K(x) for any [37], with the following estimated convergence speedWe now consider the arbitrary distribution function, P(x) and carrying capacity, K(x), and we have the following result.

Theorem 2

If the rational function , where k is a constant then is an equilibrium solution of (14)-(16). Moreover, for any the solution v(t, x) of (1)-(3) with converges to , with the estimated convergence speed as

Define the constantsfor some fixed moment . Note that since by assumption and for any by Theorem 1.

Let and be the solutions of (14)-(16) satisfying initial conditions and , respectively. According to Definition A2 the function satisfies (57) and it is an upper solution of (14)-(16). The function is a lower solution of the problem (14)-(16) with the initial condition at the moment T which is v(T, x). Then, by Theorem 1 we have and for any such thatThus, for any positive , the above inequalities implyThe following two inequalities are deduced from the same idea as above and henceCombining (18) and (19) leads toThus, integrating both sides yields.We now define and then by (19) and substituting in (14)-(16), we haveIntegrating both sides of (21) over and using the Neumann boundary conditions with the application of the Gauss theorem yieldsSince is an upper solution of (14)-(16), we haveIntegrating over the domain, using boundary conditions, and we getSince , for a real constant , it is possible to write and the gradient, will no longer be zero. Then for some and P(x), which yields the following inequalities from (22) and we haveTherefore by using Grönwall’s lemma A4, we getwhich impliesWhile we consider the lower solution, , similarly, it is easy to establish the following integrals relation.Finally from (20), (24), and (25) it follows thatwhereIf t is too large and tends to thenfor any non-negative and not identically zero initial value in . Equality (26) implies the uniqueness of the positive equilibrium solution.

For rational functions, all non-negative and non-trivial solutions converge to the positive steady-state K(x), and we can estimate the convergence rate. The result is prescribing in the following Theorem.

Theorem 3

If the function in (14)-(16) then for any the solution v(t, x) of (1)-(3) with converges to K(x), with the convergence speed estimated as

As far as the function K(x) and P(x) are rational, the proof is very straightforward and similar to Theorem 2 and is omitted.

periodicity and unique periodic attractor

When the resource function K(t, x) and positive distribution P(t, x) are both space and time-dependent, generally, there is no more a positive equilibrium solution of (1)-(3) without having the effect of harvesting. However, if P(t, x), K(t, x) and r(t, x) are periodic bounded functions with the same period , under certain conditions, we can establish the existence of a positive attracting periodic solution.

Recall again the substitution , where P is assumed to be strictly positive, bounded, Hölder continuous in x, periodic and at least in t. Then the problem (1)-(3) with becomeswhere .

We also denote and for any bounded continuous function g(t, x).

Theorem 4

Suppose that andwhere H is defined in (7), then there exists a unique periodic solution of (27)-(29) in the interval

DenoteIt is easy to check thatThese inequalities imply that and are the constant lower and upper solutions of (27)-(29) with the initial conditions and , respectively. Then according to Lemma A4 from Appendix there exists a pair of -periodic solutions and satisfyingNow we will show that if condition (30) holds. Subtracting (32) from (33), we constructIntroducing in the above equation, we haveRearrangingMultiplying both sides of (37) by and rearranging, we haveDue to periodicity of , and , the function s(t, x) is also -periodic, thereforeIn the first inequality of the previous line we used ; the last inequality follows from the condition (30). Therefore and we obtain a unique periodic solution of (27)-(29) in the interval . Moreover, according to Lemma A4 for any initial function the corresponding solution of (27)-(29) satisfies as for .

Persistence and extinction with harvesting

We now consider the harvesting term effect (i.e., ) in the system (1)-(3) to explore some advanced study connected between intrinsic growth rate, carrying capacity, and harvesting function.

Introducing in system (1)-(3) becomesThe system (39)-(41) also has a unique positive solution, and we have the following result.

Theorem 5

Let in , then there exists a unique positive solution v(t, x) of the problem (1)-(3).

The proof of the theorem follows using after the immediate substitution of (1)-(3), which brings out (39)-(41) and as a consequence of the proof of Theorem 1. According to Theorem 1, to prove it, we only need to construct an upper-lower solution pair. Consider a constant C as in (9) then which can be found according to (h2). Furthermore, choose and let which is valid since is bounded in and K, and P are bounded from below. ThenFurther,Therefore, is an upper solution of (39)-(41) by Definition A2. The function is obviously a lower solution. The right-hand side function of (39)-(41)is continuously differentiable with respect to w and we can denote the maximal derivative of g in w for each (t, x) such thatThen, the Lipschitz condition (59) holds and by Lemma A2 there exists a unique solution of the problem (1)-(3) satisfying . The positivity of the solution is obvious according to Theorem 1.

In the following stairs, we will consider the case where all the functions are time-independent, and we assume , , , and . Thus, the system (1)-(3) becomesTo study the problem (43)-(45), we first consider the following two propositions for a general parabolic equation as described in [17].

Proposition 1

( [17], Proposition 3.1) Consider the following general form of a modelAssume that for , where is a bounded measurable function if is a gradient, and if is not a gradient. If the largest eigenvalue ofis negative, then (46)-(47) has no positive equilibria, and all non-negative solutions decay exponentially to zero as .

Proposition 2

([17], Proposition 3.2, and 3.3) Suppose that with g(x, v) of class in v and in x and there exists a such that for . If the principal eigenvalue is positive in the problemthen (46)-(47) has a minimal positive equilibrium , and all solutions to (46)-(47) which are initially positive on an open subset of are eventually bounded below by orbits which increase toward as . Moreover, if g(x, v) is strictly decreasing in v for then the minimal positive equilibrium is the only positive equilibrium for (46)-(47).

Now we will demonstrate sufficient conditions for extinction and persistence for the harvesting model (43)-(45).

Theorem 6

Let the function satisfies the assumptions (h1)-(h2) and let for all . Then for any the solution of system (43)-(45) converges to zero as .

Once again, we assume and then Eq. (43)-(45) becomesRewrite the term as:Therefore, comparing (50)-(52) with (46)-(47), it shows that , , , and .

DenoteThen the linearized eigenvalue problem of (50)-(52) around the solution w(t, x) is given byTo prove the result, it is enough to show that the principal eigenvalue of (53)-(54) is negative. Multiplying both sides of (53) by P(x), and integrating over using the boundary condition, sequentially, we getTherefore, the principal eigenvalue isFrom the initial assumption, we have such thatHence, , which concludes the proof.

Theorem 7

Let the function hold the conditions (h1)-(h3) and let for . This guarantees a unique positive equilibrium state of (43)-(45) such that , , as t approaches to the infinity, v(t, x) converges to the .

If we compare the system (50)-(52) which is equivalent to (43)-(45), with system (46)-(47), then we get the linearized system around the solution w(t, x) given bywhere,By integrating and using the boundary condition we getWhich gives us the principal eigenvalueSince g(x, w) is strictly positive from the initial assumption that for all , we can conclude that .

Also,using the properties (h2) there exists K such thatHence the function g(x, w) is decreasing in w for . By comparing with Proposition 2, there exists a unique positive equilibrium solution of (50)-(52) as equivalence to of (43)-(45).

To show that, the equilibrium solution is bounded, we can use upper-lower solution pair according to Lemma A2 as defined in Appendix A. The proof is very similar to Theorem 5. For simplicity we will use model (50)-(52). To construct an upper solution, consider a constant c such thatthen . Further, choose and let which is valid since is bounded in and K and P are bounded below in . ThenAlso,Therefore is an upper solution of (50)-(52) by Definition A2. The function is obviously a lower solution. Since the steady-state is unique and bounded by upper and lower solution, so it is asymptotically stable, i.e., as for any .

Numerical simulations

In this section we perform several numerical experiments with and without harvesting effect in the model. For numerical computation, we discretize the Eqs. (1)-(3) using Crank-Nicolson scheme, which is unconditionally stable.

Numerical illustrations without harvesting

Example 1

Assume that on with the diffusion coefficient , and the reaction coefficient . If where k is a constant, for example, , and then there is a unique solution of (1)-(3).

If we compare the uniqueness for various initial conditions and different growth rates, the unique solution of (1)-(3) is valid as designed in Fig. 1. This figure justified the result of Theorem 1 for proportional functions P(t, x) and K(t, x).

Fig. 1

Average solutions of (1)-(3) using different initial values for (a) , (b) , and (c) , where , , , and

In the following example, we consider the arbitrary time-dependent function to check the positivity and uniqueness of solution for the system (1)-(3) by changing the initial densities.

Example 2

Suppose on with the diffusion coefficient , and reaction coefficient . If , there is a unique solution of (1)-(3) as shown in Fig. 2, where we choose various positive initial values to justify the uniqueness of solutions.

Fig. 2

Solutions of (1)-(3) a over the domain, and b over time, for , , , , , and for the varying initial conditions

Graphical structure of these diagrams (Fig. 2) verified the Theorem 1 for random, non-proportional and time-periodic functions P(t, x), and K(t, x).

In both Examples 1 and 2 (Figs. 1 and 2), Theorem 1 is verified either for proportional or non-proportional relations of the distribution function P(t, x) and the carrying capacity K(t, x), i.e., the bounded H defined in (7) does not affect the main result of the Theorem 1.

For Example 3, we consider the spatially distributed rational function to show the effects of different parameters.

Example 3

In this example, we consider on with the diffusion coefficient , and reaction coefficient . If , and there is a unique solution of (14)-(16) as designed in Fig. 3 due to the outcome of diffusion coefficients and multiple initial values.

Fig. 3

Average solutions of (14)-(16) over time for , , , and in a convergence study for , and , b effects of diffusion coefficients where , and c for different initial values with

Figure 3 ensures the result as justified in Theorem 3 for spatially distributed rational functions P(x) and K(x).

Since there exists a unique solution that converges to K(x), which is also true if we vary diffusion coefficient, d and initial density, (Fig. 3). Also, the conclusion of the Theorem 3 is stable for varying reaction coefficients, as depicted in Fig. 4.

Fig. 4

Solutions of (14)-(16) (a) over the domain, and (b) average solutions with respect to time for , , , and with the variation of reaction coefficients , and .

Finally, from those Figs. (3 and 4), we see that for proportionally changes distribution function if we alter diffusion coefficient d, initial density , and reaction coefficient r(x), the behavior of solution does not show any significant change and always converges to K(x) at .

At this stage, our interest is to consider the arbitrary function to observe and check the validity of theoretical results by changing the parametric values.

Example 4

In this problem, we consider on with the diffusion coefficient , reaction coefficient , and . We assume the carrying capacity , and . There is a unique solution of (14)-(16), which justifies Theorem 2 for the non-rational functions P(x) and K(x) (Fig. 5).

Fig. 5

Solutions of (14)-(16) a over the domain and b average solutions with time for , , , , and

Since it is not necessary here that the solution converges to K(x), it might converge to K(x) or not depending on the environment, which is also true if we vary the diffusion coefficient (Fig. 6(a)); specific growth rate (Fig. 6(b)) and the initial density (Fig. 6(c)).

Fig. 6

Average solutions of (14)-(16) over time for , , , and in a effects of diffusion coefficients where , b for different growth rates when , and c for different initial values with

Observing all figures (Fig. 6), it concludes that for arbitrary distribution function, if we vary the initial density of , then the solution doesn’t diverge. If we choose a different diffusion coefficient d, then a slight change will be noticed. Lastly, if we dealt with the variation of reaction coefficient r(x), the behavior of the solution changes rapidly, showing the speed of convergence.

Numerical illustrations with harvesting

Example 5

In model (43)-(45), we consider the spatial functions , , the diffusion coefficient , and initial values . The result in Fig. 7 follows the analytical behavior. We consider in (a) with various constant and non-constant harvesting functions, , and in (b) with , and .

Fig. 7

Average solutions of (43)-(45) over time for , , and with

When we allow harvesting, as in (43)-(45), numerical study shows that if the growth rate per capita is lower than the harvesting, either spatially distributed or constant, the density of average population levels eventually dropped to zero as time grows. On the other scenarios, if the harvesting coefficient is lower than the specific growth rate, the population has a stable state and this fact is analytically proved in Theorems 6 and 7.

Example 6

Similarly, in (1)-(3), we consider the time-varying rational function, ; where , growth rate , diffusion coefficient , and initial density . In Fig. 8, we consider harvesting coefficient in (a) ; in (b) ; and in (c) . This example also shows that the population is in extinction if the harvesting rate is high, see Figure 8(a), whereas the population persists periodically for higher growth rate, see for examples in Fig. 8b and 8c.

Fig. 8

Average solutions of (1)-(3) over time for , and on the domain , , and with for a ; b ; c

Example 7

In (1)-(3), we now assume the time-periodic non-rational parameters such that , , , with constant diffusion rate , and initial density . It is seen that the graphs follow the theoretical results, see Fig. 9 for (a) , (b) , and (c) . By choosing all functions in both space-time periodic, we conclude that there exists a positive solution that has been proven theoretically in Theorem 5.

Fig. 9

Periodic solution of (1)-(3) with time for , , , , , and initial value with a , b , and c

Conclusion

In this paper, we studied a reaction-diffusion model for single species with logistic growth function, where carrying capacity, K and distribution function, P are time dependent. We showed that for the described model, there exists unique solution (whether the distribution function P(t, x) is proportional to carrying capacity K(t, x) or not) and it is positive by using upper and lower solution method. The model is also valid for time independent carrying capacity. If the distribution function P(x) is proportional to K(x) then there exists unique solution v(t, x) which converges to K(x) with the convergence speed . On the other hand if K(x) and P(x) are non-proportional then it is not guaranteed that the solution will converge to K(x). If the carrying capacity, K(t, x) and distribution function, P(t, x) is a time-periodic function, then we showed that there has a unique periodic solution in some intervals with a threshhold dynamics. At the inclusion of the harvesting term, we studied the modified model and established several results connected with the time-spatial specific growth rate and harvesting function. The results reflect the natural phenomena of various species, as introduced in the earlier section of this study. Some numerical results were illustrated to justify the analytic results in a non-empty open domain for various parametric values with different initial conditions.

In the COVID era, with the initiative of One Health Research Coalition [43], the importance of the transdiciplinary research across environmental sciences and public health is highly encouraged. Recent reports show that COVID-19 can be spread to animals from human, and vice-versa [44]. In the future work, we plan to incorporate the COVID-19 models to our model following the recent works [45-49] to have the pandemic impact on the ecological results. To model and analyze the density of the species using fractional order derivative following the recent investigations in [50-53] and along with the idea herein could be a next research avenue.