Abstract

Two time-varying Beverton–Holt models are investigated in which the population of the same species evolves jointly in two coupled habitats which can be subject to population exchanges. Both habitats can have different parameterizations concerning their intrinsic growth rates and their environment carrying capacities due to different environmental conditions. Mutual fluxes of populations in-between both habitats are possible together with harvesting actions. In one of the models harvesting acts on juvenile individuals. In the other proposed model, harvesting takes place on the adult populations after the reproduction cycle they have performed has ended. The second investigated model, contrarily to the first one, relies on an “a posteriori” harvesting action to the reproductive stage which is able to modify the stocks of population. The considered harvesting can also be negative to describe repopulation actions. The equilibrium points in steady-state and their stability properties as well as the extinction conditions and the boundedness, oscillation issues, and positivity of the solutions are also investigated.

1. Introduction

Discrete-time Beverton–Holt equation-based models, based on the logistic equation, are a common tool for modelling the evolution of species which reproduce by eggs such as insects, birds, and fishes [14]. The most elementary Beverton–Holt equation is parameterized by the carrying capacity of the environment and the species intrinsic growth rate. The above parameters can be generalized to sequences in order to consider different conditions on the population evolution in different seasonal periods. Additional parameters such as the independent consumption, associated with eventual disturbances, and harvesting, associated with fishing or hunting, are also included in generalized versions of the equation. The elementary invariant Beverton–Holt equation possesses two equilibrium points, namely, the extinction one and the carrying capacity level. Under typical evolution conditions of intrinsic growth rates exceeding unity, the extinction equilibrium point is unstable while the nonextinction one is globally asymptotically stable [4, 5], implying also the population persistence. If the Beverton–Holt equation is subject to periodic carrying capacities and intrinsic growth rates of the same period then the averaged periodic sequence of population numbers lies below the average of the carrying capacities. The above, so-called, Cushing–Henson conjecture has been rigorously proved to be true by Stevic [2]. Some of their extensions for the q-difference Beverton–Holt equation have been described and proved in [3]. On the other hand, control theory issues for Beverton–Holt have been investigated and discussed in [46]. In particular, control actions have been proposed consisting of an online design of the environment carrying capacity for monitoring its suitable sequence of values necessary to match a prescribed suitable population evolution sequence. This method is claimed to be potentially applicable in semiopen environments such as some fisheries. It can be pointed out that the intrinsic growth rate and the environment carrying capacity can be related to each other, in real situations, as discussed in [7]. On the other hand, it has also been proposed an impulsive extended competition Beverton–Holt model among species from the stability point of view [8]. Beverton–Holt models as well as other mathematical related models like, for instance, Ricker’s-type models are often used in Maritime Biology for population evolution studies and population management. See, for instance [9]. In [68, 1013], harvesting actions have been investigated in the framework of extended Beverton–Holt models. Typically, harvesting is related to legal regulated fishing or hunting but illegal poaching might also be included [6, 8]. Periodicity of the solutions of the Beverton–Holt equation due to seasonality has been focused on in [14, 15]. On the other hand, a global dynamics analysis of extended models [16], as well as the appearance of bifurcations [17, 18] or the presence of resonances [19], has also been investigated. More recently, an extended Beverton–Holt model defined on isolated time scales has been discussed in [20] while an extended Beverton–Holt equation that considers discrete delays has been described in [21]. Also, more sophisticated Beverton–Holt-based models are often utilized for monitoring fishing stocks and migration flows in order to elucidate the recommended top of captures [22, 23]. The so-called independent consumption, that is, the unforeseen changes of populations due to perturbations can be characterized also under stochastic formulations. See, for instance [24, 25] and some of the references therein. Recently, a dynamic version of the discrete Ricker and Beverton–Holt model has been proposed and discussed in [26]. The population evolves according to a density-dependent population contribution and limited replenishment of the stock. The persistence of the population and the stability properties has been also investigated. Also, it is known that the exchanges of populations are useful for recovery actions of the habitats [27], what is of interest for its regeneration and its eventual exploitation.

On the other hand, evolutionary frameworks for the description and mathematical characterization of the species evolution have received important attention from long time ago. The pioneering work later on leading to those frameworks in more recent times was undoubtedly Darwin´s book “On the origin of the species” [28], which describes the fundamental evolution principle based on three axioms of natural selection related to variation, inheritance and competition. See also, for instance [29, 30]. Variation refers to different traits of phenotype of individuals belonging to the same population. Roughly speaking, inheritance is the mixture of both parents traits transmitted to the new-borns. Competition primarily refers to the fact that the best transmitted genetic inheritance to the environment make that individuals endowed with them can better survive related to the others. See, for instance [2831], and some related references there in. Competition can also refer to the case of various species sharing a certain habitat which are competing for food, refuge for defence of predators, etc. and to the case of predator-prey competition for survival. It can be pointed out that, in the background literature, it has been also usual to develop Ricker’s-type models for the characterization of population evolution in the frameworks of interspecific and intraspecific competition related, respectively, to competition between individuals or groups of the same or different species. See, for instance [32, 33], although the related literature is very rich.

Also, Beverton–Holt-type evolution models are also useful for those purposes. In the above context, a discrete evolutionary Beverton–Holt model has been proposed and studied in [29], supported by the evolutionary game framework by considering the strong Allee effect related to predation saturation. Also, Neimar–Sacker bifurcation with chaotic behaviour has been detected in the model of this research, while its asymptotic stability is studied from fixed point theory. Also, a new approach for a specific class of discrete time evolution models has been proposed in [30] for a class of new single- and multispecies evolutionary competition models by considering the model as a nonautonomous difference equation.

In this paper, two time-varying Beverton–Holt models are proposed and studied which consider that individuals of the same species are living and reproducing in two different habitats which are not mutually isolated since they admit populations interchanges. Those habitats can have different characteristics of intrinsic growth rates and environment carrying capacities due, for instance, to be subject to disposal, different characteristics of temperature, humidity, available refuge, etc. Mutual flows of population in-between both habitats are considered as well as the existence of harvesting actions, basically, legal hunting or fishing, and/or illegal poaching, and also eventual repopulation actions, those ones being formulated as negative harvesting. The migrations between habitats as well as the harvesting actions are considered to be proportional to the numbers of individuals while the whole balances of changes of population they generate with respect to the standard population variations, being governed by the intrinsic growth rates and the carrying capacities, are considered as a generalized harvesting contribution. It can be pointed out that exchanges of populations could be useful, in particular, in the context, for instance, of habitat recovery actions, which are receiving certain interest in part due to the effects of climate changes. See, for instance [27], and some of the references therein. The first proposed model considers that the generalized harvesting takes place on juvenile individuals so that they do not contribute to the standard population growth at the reproductive stage. In typical real cases, this problem view could describe, for instance, repopulations with either fingerlings or juvenile incorporations to the habitats or by removal of some of them for later exploitations in fisheries. In this sense, the generalized harvesting is viewed as an “a priori” action on the nonadult population which can coexist with the reproductive stage performed by adults. The second model considers that the generalized harvesting takes place on the adult populations after the reproduction cycle they have performed has ended. It can be viewed as an “a posteriori” generalized harvesting action to the reproductive stage which is able to modify the stocks of population.

This paper presents a contribution to the scientific literature by introducing two interconnected time-varying Beverton–Holt models with the inclusion of generalized harvesting actions alongside classical reproduction stages. Through a rigorous analysis, we establish the well-posedness of both models, thoroughly explore the stability of equilibrium points, and delve into the oscillation properties of the solutions. Beyond its theoretical significance, this research holds practical promise for the management of production fisheries operating under dynamic and diverse conditions, offering valuable insights that can inform more effective decision-making in this critical domain.

The paper is organized as follows. Section 2 establishes both models based on two time-varying Beverton–Holt coupled equations being, respectively, eventually subject to generalized “a priori” and “a posteriori” generalized harvesting actions, that is, in parallel to the standard reproduction stage or after such a stage has ended. Section 2 also gives sufficiency-type nonnegativity conditions for the solution of both models as well as the characterization of the equilibrium points of both models in stationary regime. The local stability and the oscillation conditions of Model 2 are also investigated. On the other hand, the characterization and the local stability properties of the equilibrium points and some oscillation properties of Model 1 are investigated in Section 3. The equilibrium points are shown to be identical in both models although their stability conditions differ. The boundedness of the solutions, under reasonable conditions on the time-varying sequences of parameters which parameterize the models, is also proved in Sections 2 and 3. This fact concludes that, under such conditions, global stability is also achieved. Some numerical simulations are given and discussed in Section 4, and, finally, conclusions end the paper.

1.1. Notation

and denote, respectively, the sets of nonnegative and nonpositive integer numbers.

and denote, respectively, the sets of positive and negative integer numbers.

and denote, respectively, the sets of nonnegative and nonpositive real numbers.

and denote, respectively, the sets of positive and negative real numbers.

and stand, respectively, for a matrix with at least one positive entry, respectively all positive entries.

2. Problem Statement and Main Results

The parameterizations of Beverton–Holt equations can be given by sequences of intrinsic growth rates and environment carrying capacities in a more general context that periodic sequences associated, for instance, with seasonality issues [26]. In the sequel, two time-varying discrete Beverton–Holt models are proposed and studied which consider that numbers of members of the same species are living and reproducing in two different habitats or environments. Those habitats can have different characteristics of intrinsic growth rates and environment carrying capacities due, for instance, to different characteristics of temperature, humidity, available refuge, etc. Mutual flows of populations in-between both mentioned habitats are considered as well as the existence of harvesting actions, basically, legal hunting or fishing or illegal poaching, and also eventual repopulation actions being formulated as negative harvesting. It includes also eventual independent consumption effects on the populations. The migrations between habitats as well as the harvesting actions are considered proportional to the numbers of individuals and the whole balances of changes of population they generate with respect to the standard population variations are considered as a generalized harvesting. The first model considers that the generalized harvesting takes place on juvenile individuals so that they do not contribute to the standard population growth at the evolution stage where the standard reproduction takes place. In real cases, this could describe, for instance, repopulations with fingerlings or juvenile incorporations to the habitats or removal of some of them for later exploitations in fisheries. The second model considers that the generalized harvesting takes place on adults after the reproduction cycle has ended.

2.1. Model 1 Subject to “A Priori” Generalized Harvesting on Juvenile Individuals

The population of the two groups of the same species evolves according to the following time-varying coupled Beverton–Holt equations:

under initial conditions for . The populations in the habitats at the -th sampling instant are, respectively, and . The sequences and are the intrinsic growth rates and the environment carrying capacities, respectively. The terms involve eventual couplings between both populations and describe the generalized harvesting terms which can be positive, negative, or null. The coefficients describe the fractions of the respective populations which leave each of the habitats, respectively, 1, 2 to migrate to the other habitat, respectively, . The coefficients represent the fractions of immigrated individuals to each habitat from the other habitat, some of them could not reach their final new habitat so that the constraints and hold. The terms ; , with positive coefficients, describe the standard harvesting (i.e., fishing/hunting) including the legally regulated one and the eventual illegal poaching. If a coefficient for some , and since the generalized harvesting are allocated with minus signs in (1) and (2), that means that the eventual repopulation dominates the eventual standard harvesting action of fishing or hunting.

Note that it is considered, for the sake of simplicity, that the eventual reproduction takes only place over of the existing numbers of individuals in the habitat. In this way, if there is no population in one of the habitats, the repopulation strategy is not activated and, if it is activated, then the numbers of new individuals do not exceed the already existing population. Thus, and, if , then , . Also, note from (1)–(4) that the populations are nonnegative for any nonnegative initial conditions if the generalized harvesting sequences are nonpositive. This can happen as related to repopulation. However, if the generalized harvesting sequences have positive members, those should be upper-bounded to keep the nonnegativity of the populations. The following assumption is made in the sequel related to this concern.

Assumption 1 (Limited positive harvesting sequence). , , , and are bounded sequences, andwhich can be expressed equivalently asNote that Assumption 1 is equivalent to.equivalently,The following result addresses the nonnegativity of the populations under nonnegative initial conditions. It is used that, for all and , it holds that , so that and .

Proposition 2 (Nonnegativity of the solution of Model 1). Both populations of Model 1 are nonnegative for all samples under nonnegative initial conditions if for any and eitherunder the necessary conditionwhich always holds, irrespective of , if , or holds under the necessary condition ;

Proof. From (3) and (4) in (5):equivalently,which yield directly the proof.
The total generalized harvesting in both habitats is given bysince and .Thus, implies that so that at least one of the is negative implying a dominating repopulation action over alternative population interchanges or fishing/hunting in the corresponding habitat. The inverse populations evolve according towhereNote that, for , as expected. Note also that, as expected as well, if then so that a positive (respectively, negative) generalized harvesting leads to a negative (respectively, positive) contribution on the corresponding inverse population dynamics for any since, if and , one has even if thatwhile the respective contributions of and to and are defined with opposed signs in the corresponding formulas.
A further assumption for the Beverton–Holt equation is the following one which is inspired in a usual parallel one for the case of one habitat for the species. For that single habitat case, an intrinsic growth rate exceeding the unity threshold guarantees the species nonextinction.

Assumption 3 (Intrinsic growth rate exceeding the unity threshold). Assumption 1 holds with ; , .
In summary, the standard nonnegativity conditions for Model 1 are , together with constraint (5) to be fulfilled by the generalized harvesting sequence.
Note that by weakening the above conditions, one gets that if and then if . One also concludes that the -th population is zero, implying extinction, if . Note that in the proposed model the carrying capacities and intrinsic growth rates can be, in general, distinct for both habitats.

Remark 4. Note that the generalized harvesting contributions to the -th sample from the -the sample in Model 1 assume implicitly that the involved individuals in hunting or fishing of migration or immigration do not contribute to reproduction at the -th sample. Typical reasons are that they are moving due to migration, they are being eliminated from the population due to hunting or fishing and/or they are introduced being juvenile individuals for stock repopulation but they are not still valid for the species reproduction. Note also that the model relies on the fact that the sampling instants are not instantaneous while they have a certain potentially nonnull duration , that is, the sampling instants are ; while the intersample period fulfils .

2.2. Model 2 Subject to “A Posteriori” Generalized Harvesting on Adult Individuals

Compared to Model 1, it is now assumed that the generalized harvesting occurs at the -th sample. That is fishing/hunting and/or migrations which can alter the nominal recruitment (i.e., that evolving for null harvesting) occur after the species reproduction period. The harvesting function now, depends on the -th sample in order to accurately depict the harvesting occurring on adult individuals after the reproduction stage. Thus, if the Beverton–Holt equation represents the effect of the reproduction at the k-the sample on the at the -th sample, and since harvesting occurs thereafter, this sequence must be defined at -th sample since it does not affect to the reproduction stage at the -th sample. Contrarily in Model 1, harvesting has been assumed on juvenile individuals so that it would affect to the reproduction stage at the next -th. Therefore, harvesting has been scheduled at the -th sample for Model 1. Thus, the proposed Model 2 becomes: under initial conditions for . The above equations can be rearranged as follows:where

Remark 5. It is assumed in Model 1 that the harvesting actions are performed on juvenile elements which migrate of their habitat or either they are introduced, for instance, via repopulation, or removed, for instance, under controlled recovery of juvenile elements for later breeding in fisheries or for selling for consumption as it happens with fishing actions on young eels. So, these numbers do not contribute actively to reproduction through the parameterization of the intrinsic growth rate. In Model 2, generalized harvesting, such as fishing/hunting, repopulation, or migration actions, take place on adults after the reproduction period took place. For presentation convenience, it is first discussed Model 2 along the next section.

2.3. Nonnegativity of the Solution of Model 2

Assumption 1 or its stronger version Assumption 3 given on Model 1 are also kept for Model 2 while in addition, one assumes for Model 2 that:

Assumption 6 (Moderate negative harvesting). ; .
Assumption 6 implies that negative harvesting, for instance, repopulation (of adults) after the reproduction period at some sampling instants requires for its modulus to be less than ; .
The upper-bounding constraint (5) on the generalized harvesting for nonnegativity of the solution now takes the form.The nonnegativity of the solution of Model 2 for finite initial nonnegative conditions is characterized in the subsequent result. It is also discussed the conditions from previous samples to the current to be nonnegative provided that the solution at the current sample is nonnegative.

Proposition 7 (Nonnegativity of the solution of Model 2). The following properties hold:(i)Assume that, for each , any of the two following sets of constraints holds:(a)Migrations between both habitats with harvesting actions consisting of repopulations:(b)Null immigrations between both habitats with either repopulation or moderate fishing/hunting:Then, imply that ; .(ii)Assume that for each the immigration factors satisfy the constraints:Then, for some arbitrary implies that ; .

Proof. First, note that the coefficient matrix of the left-hand-side of (21) is nonsingular if and only ifThen,wherewherethenNote that for each , , , and is positive harvesting (i.e., fishing or hunting) while indicates repopulation. If ; then a part of the migrated individuals from each of the habi1tat does not reach successfully the other habitat. Then,(a)Conditions (25) imply that, for any , imply that with repopulations () and nonnull migration/immigration. In this case, , i.e., the adjoint matrix has strictly positive entries) and leading to , that is, is strictly positive.(b)Conditions (26) imply that, for any , imply that . In this case, , i.e., the adjoint matrix has nonnegative entries with the diagonal ones being positive) and leading to .Property (i) has been proved. Property (ii) follows since one has that is nonsingular andThis implies that for any and , for if for any and .
Note that the constraint of Property (ii) of Proposition 7 implies that

2.4. Boundedness of the Solution of Model 2

It is now discussed the boundedness of the solution of Model 2 under some weak reasonable constraints:

Proposition 8 (Boundedness of the solution). Assume that the parametrical sequences , , and are bounded for and that for . Assume that Assumption 1, with eventually its stronger version Assumptions 3 and 6 hold. Assume, furthermore, that both habitats are mutually coupled in the sense that for infinitely many samples. Then, the sequences and are bounded for any given finite nonnegative initial conditions.

Proof. Proceed by contradiction arguments. Note from (22) thatwith .Thus, if either as and for , or if for , then one concludes in both cases that for and also that as . Therefore, and for is not feasible so that if one of the populations diverges the other one diverges too so that their solutions are jointly unbounded sequences. Now, it follows from (32) that the matrix sequence as .Then, the spectral radius of converges to zero as ,from (32), while it is less than for any prefixed real constant and all nonnegative integer number for some finite . Therefore, if is finite then, from (29),For any given finite initial conditions. Hence, one reaches a contradiction to as . Therefore, .

2.5. Equilibrium Points of Model 2

Next, the eventual equilibrium points of Model 2 are characterized in the case when the parameterization sequences converge to constant values. It will be seen that the discussion of the feasibility of the nonextinction equilibrium points, in terms of their positivity, is not a trivial task since, in the general case of existence of migration and harvesting, they are the nonunique solutions of a cubic algebraic equation. By introducing some reasonable constraints on the limits of the parameterizing sequences, the existence of at least a feasible nonextinction equilibrium point is proved, and this fact is then corroborated in the numerical simulations described in Section 4.

Assume that all the parametrical sequences in Model 2 are constant and denote the eventual components of the equilibrium point as and which are the equilibrium points of each of the two habitats. Note the following features:(a)From (17)–(20), is the joint extinction point of both populations which implies also that null generalized harvestings as well as null population interchanges due to mutual migrations as well as null fishing/hunting actions.(b)From (17)–(20), Thus,(b1) if then the extinction of implies that of and, if , then the extinction of implies that of (b2) if (that is, if there are joint mutual migrations from each habitat to the other habitat) then iff (b3) if then the extinction of does not necessary imply that of , but the extinction of can happen, and, if , then the extinction of does not necessarily imply that of although such an extinction can happen(c)The combination of (17)–(20) for a stationary solution yields the following constraints:

Now, if then if and, equivalently if . Then, if both and are nonnull then.with so that implies that(1)Either and . Since , the second constraint is further constrained to ;(2)Or and which is impossible for , that is if, .

A similar discussion can be applied for . The case of numerator and denominator of the second expression in (40) being negative is not feasible since it leads to the contradiction if

If both numerator and denominator of that expression are positive then under the constraint:

Note that if then for . This means that the equilibrium points in both habitats equalize the respective environment carrying capacities if the mutual interchanges of population coming from the other habitat equalize. Equation (38) may be rewritten as

Since and for the nonextinction equilibria then.

By zeroing the first member of the above relations, one obtainsubject to leading to the necessary condition:

By zeroing the second member of the above relations, one obtainssubject to leading to the necessary condition:

The given derivations are compacted as follows:

Proposition 9. The jointly nonextinction stationary values for both species satisfy the constraint:

Another related result for the steady states follows:

Proposition 10 (Constraints for nonzero equilibrium points). The following properties hold:(i)The jointly nonextinction stationary values for both species satisfy the constraint:(ii)If , then,subject to , andIf , then,subject to , and(iii)The components of the nonzero, and, in general nonunique, nonextinction equilibrium point satisfy the constraints:with , ; with and for some .

Proof. The equilibrium conditions might be expressed asA nonzero equilibrium vector is subject to the condition that the above coefficient matrix be rank defective, that is, . Property (i) has been proved.
To prove Property (ii), note that, if then so that subject to . Furthermore, and then it has to be solved in for given the following equation:whose positive root is the nonzero equilibrium point of the second habitat if .
Similarly, if then so that subject to and proceeding in a similar way as above, we get the equilibrium point component if as the positive zero of the second order equation in for the given equilibrium component :which completes the proof of Property (ii).
To prove Property (iii), note that the extinction point satisfies (50) and also nonzero equilibrium points have to satisfy (50) so that for some nonnegative scalars and with , the subsequent constraints have to hold in order for the null space of the coefficient matrix not to be identically zero:which lead towith.
, ; and with for some real . However, note that if or while then the solution sequence would not be unique so that (61) has to be fulfilled for calculated from (61) with ; which leads to (55).
Some more general results follow below without considering the constraint that of Proposition 10.

Proposition 11. Let and . Then, for any nonzero real , one has:

Proof. Note by equalizing to zero the two first amounts of (38) that the following constraints holds for nonzero and :andThus, in order for (64)–(65) to have the same solutions in , the coefficients of and the independent term have to be mutually proportional with the same nonzero proportionality constant satisfying.Then, for any nonzero real , one has:leading to (62)–(63).

Proposition 12 (Nonextinction consensus). The equilibrium point components are identical and nonzero ifsubject to ; and it holds

Proof. On gets from (38) that if thenwhich is equivalent towhich is compacted as (68) subject to ; and (69).
The local stability of the equilibrium points is now discussed via the convergence properties of the Jacobian matrix at such points.

Proposition 13 (Local stability of the equilibrium points). The equilibrium point of Model 2 is locally asymptotically stable ifequivalently, if the intrinsic growth rates are lower-bounded, related to the remaining stationary parameters, accordingly to:

In particular, the extinction equilibrium point is locally asymptotically stable if the intrinsic growth rates are upper-bounded, related to the remaining stationary parameters, accordingly to

Proof. Consider (21)–(23), which lead to (32), with ; ; , for a constant parameterization at the asymptotic limits of the parameterizing sequences , , , , and for ; . The limit evolution equation (22) for the limits of the parametrizing sequences can be expressed asThe Jacobian matrix at the equilibrium point becomeswhich leads toafter replacing the values of the equilibrium points from (55). Since the spectral radius of a square matrix is less than or equal to any of its norms then, by taking the and the norms in (77), the spectral radius of the Jacobian satisfiesequivalently, if (72) or its equivalent (73) hold. The conditions (74)–(75) are got as particular cases of local asymptotic stability if since in that case the Jacobian matrix becomes to beso thatguaranteed by (74).

Remark 14 (Oscillation condition). The oscillation condition for a period of samples is easy to characterize from (32). Note that if there exists a set of values of the solution ; ; for a parameterization , , , , ; ; , then such values define an oscillation of the solution of period if for , whereand is the second-order identity matrix.

3. Further Stability and Oscillation Properties of Model 1

The equilibrium points of Model 1 are the same as those of Model 2 since they satisfy also the same stationary constraints , with , , , , and for ; . Then, both generalized harvesting sequences take identical steady-state values at the equilibrium points. However, their evolution dynamics are distinct and the local asymptotic stability conditions of the equilibrium points differ as well. Therefore, there are some variations on the sufficiency-type conditions related to guaranteeing the basic properties. Note that, for Model 1, equation (32) becomes modified as follows:

Proposition 15 (Local stability of the equilibrium points). The equilibrium point of Model 1 is locally asymptotically stable if and only if the spectral radius of the Jacobian matrix at defined byis less than one while the evolution sequence of populations is nonnegative, that is, if , andor

In particular, the extinction equilibrium point is locally asymptotically stable if for , or if .

Proof. The equilibrium point is locally asymptotically stable if the spectral radius of the Jacobian matrix is less than one while the evolution sequence of populations is nonnegative. By simple inspection, both conditions hold if (84) holds or if (85) holds under similar arguments as in the proof of Proposition 13. The extinction equilibrium point is locally asymptotically stable if for or if .

Remark 16. Note that, if there is just an habitat under consideration, say the habitat 1, then the generalized harvesting does not include interhabitat migrations so that and the equilibrium points are (extinction) and if . See, for instance [46]. It turns out that is locally asymptotically stable if since while is locally asymptotically stable if since . If there is standard harvesting under the forms of fishing/hunting or repopulation (the above lower bound ensures positivity) then the nonextinction equilibrium point is locally asymptotically stable if . On the other hand, note from Propositions 13 and 15 that, if for , then the extinction equilibrium is locally asymptotically stable for both models.

Remark 17 (Oscillation condition). If there exists a set of values of the solution ; ; for a parameterization , , , , ; ; , then such values define an oscillation of the solution of period if for , where

Proposition 18 (Solution boundedness). If with for any and any such that , then for is bounded if is finite for provided that , where if and if for . A sufficient condition for that is that with for i ∈ .

Proof. Consider three cases:

Case 19. if ; , thenunder nonnegativity of the solution, where, if , if , provided that ; , . Since if ; , , thus, if as for some , then one has from (87) thatThe result has been proved if the intrinsic growth rate exceeds unity.

Case 20. Assume that . Then, for each and any either or, if , thenso that for is bounded if is finite for .

Case 21. Assume that . Then,where if and if and
for .

Remark 22 (Global stability). Note that Proposition 8, related to Model 2, and Proposition 18, related to Model 1, conclude that, under reasonable conditions of the time-varying sequences which parameterize the models, the global stability is also achievable since the solution sequences are globally bounded for all samples.

4. Numerical Simulations

This section contains three examples that illustrate some of the main theoretical results proved and discussed in the previous sections:

4.1. Example 1

This example illustrates the results of Proposition 2. The sequences , , , and , for , are, respectively, generated by means of the following difference equations: with the following values for the parameters:related to the habitat 1 and the parameters:related to the habitat 2. The following initial conditionsare considered for obtaining the evolution of the sequences generated by (91). Namely, the resulting sequences arefor and . In this way, for instance, the sequence starts with the value and it converges to the value as k tends to infinity. Moreover, the sequence is monotonically nonincreasing since . The same result can be said about all the sequences defined by (91) from the choice of the values for the parameters in such equations and the initial conditions for , , , , and , for . Figure 1 shows the evolution of the species populations , for , if the population are initially and while Figure 2 displays the evolution of and , for , where and .

One can see in Figure 2 that , for , which guarantees the nonnegativity of the solution of Model 1 as Proposition 2 establishes provided that . In summary, Figure 1 shows that , for , , accordingly to the result of Proposition 2 since the conditions and , for , and , are satisfied.

4.2. Example 2

This example illustrates the results of Proposition 7. The sequences , , , and , for , are, respectively, generated by means of the following difference equations:with the following values for the parameters:related to the habitat 1 and the parameters:related to the habitat 2. The following initial conditionsare considered for obtaining the evolution of the sequences generated by (96). Figure 3 shows the evolution of the species populations , for , if the population are initially and while Figure 4 displays the evolution of , for , where , for , and and .

One can see in Figure 4 that:(i), which implies that , for , , and(ii), which implies that , Which guarantees the nonnegativity of the solution of Model 2 as Proposition 7(a) establishes. As a consequence, the fact that , for , is shown in Figure 3 accordingly to the result of Proposition 7.

4.3. Example 3

This example illustrates the results of Proposition 10. All the parametrical sequences associated to both habitats are considered constant, namely:

The model 2 potentially has three equilibrium points which are obtained from (17)–(20) by taking into account that for at the equilibrium point. In this way, one obtains that

From (101), with and , one obtains by direct calculations that

From (101), with and , and introducing (102), one obtains that the population of the habitat 1 at the potential equilibrium point is one of the roots of the following cubic equation:with

Only the real and strictly positive roots , for , of (103), if they exist, have significance in a population model. Then, the potential equilibrium point really exits if , for , is real and positive and the population for the habitat 2 obtained from (102) by substituting , with being real and positive, in (102) is also real and positive. Such a fact is fulfilled if:

By considering the values for the parameters in (100) the roots of the polynomial (103) are the following:

Only the value satisfies the condition (105) since . Then, there is only an equilibrium point where the population in each habitat is, respectively, given by and . Such a result can be seen in Figure 5 where the evolution of the populations of the habitats 1 and 2 are displayed in Figure 5 when the initial condition is and . This result is in accordance with the result in Proposition 10(iii) related with the convergence of the model 2 to a nonzero equilibrium point.

5. Conclusions

This paper has proposed and discussed two discrete coupled time-varying Beverton–Holt which consider that individuals of the same species are living and reproducing split into two nonisolated different habitats, in fact, admitting mutual populations interchanges, standard harvesting (fishing, hunting, illegal poaching, etc.), and negative harvesting being associated to repopulation actions. Both habitats can have different parameterizing characteristics of intrinsic growth rates and environment carrying capacities due, for instance, to their different characteristics of temperature, humidity, available refuge, etc. The migrations between habitats, as well as the various harvesting actions, are considered to be proportional to the available numbers of individuals while the whole balances of changes of population they generate with respect to the standard population variations, being governed by the intrinsic growth rates and the carrying capacities, are jointly considered as a generalized harvesting contribution which can have a positive, null, or negative global balance at each sample. It can be pointed out that exchanges of populations could be useful for preservation or exploitation, in particular, in the context, for instance, of habitat recovery actions, which are receiving certain interest in the last years. The first proposed model relies on a generalized harvesting which takes place on juvenile individuals so that they do not contribute to the standard population growth at the reproductive stage. In typical real cases, this approach might describe, for instance, repopulations with either fingerlings or juvenile incorporations to the habitats or by removal of some of them for later exploitations, for instance, in fisheries. In this sense, the generalized harvesting is viewed as an “a priori” action on the nonadult population which can coexist with the reproductive stage performed by adults. The second model considers that the generalized harvesting takes place on the adult populations after each reproduction cycle they have performed has ended. It can be viewed as an “a posteriori” action to the reproductive stage which is able to modify the stocks of population. The equilibrium points of the stationary solution in the presence and absence of harvesting action have been formally characterized as well as their local asymptotic stability properties in the case of intrinsic growth rate exceeding unity and eventual execution of harvesting actions. The equilibrium points are identical for both models although their respective stability conditions can be distinct. It has been found that three nonextinction equilibrium points can exist of which, at least one, is feasible under reasonable constraints of the immigration levels related to the remaining model parameterizing parameters. Some numerical examples have also been discussed. The restriction to two habitats of the proposed models has been done for exposition clarity but the formal extension of the results to habitats is direct.

Data Availability

No underlying data were collected or produced in this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors are grateful to the Basque Government for its support through Grant Basque Government (IT1555-22), to MCIN/AEI 269.10.13039/501100011033 for Grant (PID2021-1235430B-C21), and to MCIN/AEI 269.10.13039/501100011033 for Grant (PID2021-1235430B-C22).