applied sciences Article On the Properties of a Class of Impulsive Competition Beverton–Holt Equations Manuel De la Sen 1,2, * , Asier Ibeas 3 , Santiago Alonso-Quesada 1,2 , Aitor J. Garrido 2,4 and Izaskun Garrido 4 1 Department of Electricity and Electronics, Faculty of Science and Technology, Campus of Leioa, University of the Basque Country, 48080 Bilbao, Spain; santiago.alonso@ehu.eus 2 Institute of Research and Development of Processes, Faculty of Science and Technology, Campus of Leioa, University of the Basque Country, 48080 Bilbao, Spain; aitor.garrido@ehu.eus 3 Department of Telecommunications and Systems Engineering, Universitat Autònoma de Barcelona, 08193 Barcelona, Spain; Asier.Ibeas@uab.cat 4 Department of Automatic Control and Systems Engineering, Faculty of Engineering of Bilbao, University of the Basque Country, Rafael Moreno No. 3, 48013 Bilbao, Spain; izaskun.garrido@ehu.eus * Correspondence: manuel.delasen@ehu.eus Abstract: This paper is devoted to a type of combined impulsive discrete Beverton–Holt equations in ecology when eventual discontinuities at sampling time instants are considered. Such discontinuities could be interpreted as impulses in the corresponding continuous-time logistic equations. The set of equations involve competition-type coupled dynamics among a finite set of species. It is assumed that, in general, the intrinsic growth rates and the carrying capacities are eventually distinct for the various species. The impulsive parts of the equations are parameterized by harvesting quotas and independent consumptions which are also eventually distinct for the various species and which  control the populations’ evolution. The performed study includes the existence of extinction and  non-extinction equilibrium points, the conditions of non-negativity and boundedness of the solutions Citation: De la Sen, M.; Ibeas, A.; for given finite non-negative initial conditions and the conditions of asymptotic stability without or Alonso-Quesada, S.; Garrido, A.J.; with extinction of the solutions. Garrido, I. On the Properties of a Class of Impulsive Competition Keywords: difference equations; discrete Beverton–Holt equation; impulsive equations; competition Beverton–Holt Equations. Appl. Sci. Beverton–Holt equations; equilibrium points; non-negativity; boundedness 2021, 11, 9020. https://doi.org/ 10.3390/app11199020 Academic Editor: José A.F.O. Correia 1. Introduction Received: 12 July 2021 The single-species non-autonomous discrete Beverton–Holt equation has been widely Accepted: 17 September 2021 proposed for the study of the dynamics of populations in ecology and, in particular, for Published: 28 September 2021 the characterization of the population growth dynamics. Its typical usefulness is related to the evolution of some species which reproduce by eggs. See, for instance, [1–4] and Publisher’s Note: MDPI stays neutral some of the references therein. Some usefulness of the models has been pointed out related with regard to jurisdictional claims in to the fishery industry relying on the exploitation of species such as, for instance, coho published maps and institutional affil- salmon, plaice, haddock and others [1]. The model has also been reported to be useful for iations. describing the evolution of the acorn wood pecker. Therefore, a certain effort has been devoted to researching the mathematical properties of such an equation as, for instance, the equilibrium points and their stability properties, as well as proving some classical conjectures on Beverton–Holt equations concerning the positivity of the solutions and the Copyright: © 2021 by the authors. existence of stable oscillations [5–14]. The discretization could potentially be extended to Licensee MDPI, Basel, Switzerland. time intervals, including two samples of, in general, distinct sampling periods, to separate This article is an open access article the evolution of egg/larvae from that of reproductive adults in some species. In particular, distributed under the terms and two different sampling periods can be involved when the stages of egg/larvae evolution conditions of the Creative Commons and those of adult evolution have different relevant average time intervals. For theoretical Attribution (CC BY) license (https:// multirate sampling techniques of analysis involving two or more combined sampling creativecommons.org/licenses/by/ periods, see, for instance, [15] and some of the references therein; see also [16–22]. An 4.0/). Appl. Sci. 2021, 11, 9020. https://doi.org/10.3390/app11199020 https://www.mdpi.com/journal/applsci Appl. Sci. 2021, 11, 9020 2 of 38 important effort has been devoted to rigorously prove the so-called Cushing–Henson Conjectures [8]. The standard single-species non-autonomous Beverton–Holt Equation is µ(k)K (k) x (k) x ( k + 1) = , k ∈ N0 : = N ∪ { 0 } K ( k ) + ( µ ( k ) − 1) x ( k ) where N is the set of natural numbers, for some given x (0) = x0 ≥ 0, where {µ(k)}∞ k =0 and {K (k)}∞ k =0 are, respectively, the intrinsic growth rate sequences of the species and the sequences of the carrying capacity of the habitat. The subscript k refers to each sampling instant in the evolution of the equation. In the above equation, biologists refer to the popula- tion xk as the spawning stock at the continuous time kTs , with Ts being the sampling period in-between consecutive samples [4]. An interpretation of the Beverton–Holt equation is that it is the discrete version of the logistic equation running at a sampling period Ts which is the time interval in-between two consecutive sampling instants tk = kTs and tk+1 = (k + 1) Ts ; ∀k ∈ N0 . It is well-known that the Beverton–Holt equation, sometimes referred to as the Pielou logistic difference equation, is equivalent to the discrete Verhulst logistic equation x (k ) = r (k) x (k )/( x (k) + ν(k)); ∀k ∈ N0 , under a given initial condition x (0) = x0 [13]. The time-invariant version of the Verhulst equation is the discretization of the logistic . continuous-time time-invariant equation x (t) = r x (t)(1 − x (t)/K ); ∀t ∈ R0+ = R+ ∪ {0}, x (0) = x0 with x (t) being the size of the resource population [13]. On the other hand, Stevic [12] has proved the conjectures of Beverton–Holt analytically [8], for periodically varying carrying capacities, which have been then rigorously confirmed. The first one of such conjectures claims the existence of a periodic steady-state regime in the event that the carrying capacity varies in a periodic fashion. Such a periodic steady-state is a global attractor of all the positive solutions. The second conjecture establishes that the average value of the periodically varying values of the solution is smaller than the average value of the carrying capacity over one such sampling interval. In addition, the Beverton–Holt equation has been extended by including discontinuities at the sampling time instants. Such discontinuities might be interpreted as a presence of spawning stock impulses in the corresponding continuous-time equation. The discontinuities are parameterized by two more sequences of parameters, the so-called harvesting (typically, fishing or hunting) quota sequence {(1 − b(k))}∞ ∞ k =1 and the so-called independent consumption sequence { d ( k )}k=1 , which can include migrations “into” and “from” the studied region, namely, of positive and negative values, respectively. It could also include, in the first case, re-population by adults which is not properly reflected in the non-impulsive version of the equation since it is not related to the evolution from the larva stage to the adult one. In this way, an artificial repopulation may be interpreted in this context as a migration into the region considered for the population study. On the other hand, the loss of population by dead are included in the basic Beverton–Holt equation by appropriate settings of the environment carrying capacity and intrinsic growth rate. It can be pointed out that the regulation of the fishing quota is a very relevant problem since the fishing resources are becoming increasingly scarce. As a result, the impulsive extended Beverton–Holt model is parameterized by four parameters, or sequences, namely, the environment carrying capacity, the intrinsic growth rate, the harvesting quota and the independent consumption. It has to be clarified that, since input impulses in continuous-time dynamic systems translate into bounded discontinuities in the state (see, for instance, [23–25]), it is usual to refer to their counterpart discretized versions (with bounded steps at sampling time instants) as impulsive discrete systems. See, for instance, [13,26,27]. The Beverton–Holt equation has also been studied from a control point of view. For instance, the carrying capacity is designed in [28,29] as a control variable, within certain limiting bounds, to achieve a prefixed suitable solution. That strategy is feasible in certain closed or semi-closed habitats, such as, for instance, in fisheries or certain agricultural exploitations. There are also some more recent studies available where two-stage extended Beverton–Holt equations have been described in the sense that the solution depends on its values on two consecutive previous samples. The relevant functions of two arguments (that is, the two preceding values of the sampled Appl. Sci. 2021, 11, 9020 3 of 38 solution) that define the evolution dynamics are assumed, by different researchers, to be non-decreasing in both or one of those arguments [30,31]. The study of periodic solutions has been also studied in the non-impulsive and in the impulsive Beverton–Holt equation since some of its parameterizing parameters, such as, for instance, the carrying capacity or the harvesting quota, are often subject to cyclic variations. See, for instance, [26,32,33] and some of the references therein. The dynamics associated with strong and weak Allee effects are studied in [34]; these lead, respectively, to extinction, for population stocks under a critical density, and to recovery, after small density levels are reached without the presence of a critical value. The results are applied to the study of the coexistence of two competing species subject to Allee effects and contest competition. It can be pointed out that a Beverton–Holt equation type is proposed and studied in [35] for a set with, in general, more than two competing species. It is found that the species with the best fitness outcompete the remaining ones. For the study of the involvement of two sampling periods for two stages of evolution of mosquito, see, for instance [35] and some of the references therein. In this paper, a group of impulsive discrete Beverton–Holt equations reflecting compe- tition among a finite number of species within the same habitat is considered with eventual discontinuities of the solution at the sampling instants due to harvesting or independent consumption influences. The above-mentioned set of equations involve competition-type coupled dynamics among the involved species. It is assumed that, in general, the intrin- sic growth rates might be eventually distinct for the various species and, since they can have, in general, different fitness to the habitat, the various carrying capacities might be also distinct, in general. The impulsive effects at sampling instants in the equations are parameterized by the various harvesting quotas and the independent consumptions that can be also distinct for the various competing species. The performed study includes the existence of extinction and non-extinction equilibrium points, provided that the sequences of parameters converge to appropriate values, the conditions of non-negativity and bound- edness of the solutions, for any given sets of finite non-negative initial conditions, and the description of sufficiency-type conditions, for asymptotic stability to any equilibrium point and, in particular, for either non-extinction or extinction situations. The case of intrinsic growth rates being less than unity and implying asymptotic extinction under supplementary conditions involving the independent consumption to the harvesting quota is also discussed through the second Lyapunov stability method. The paper is organized as follows: Section 2 defines the whole impulsive competi- tion Beverton–Holt model, the equilibrium points and the positivity and boundedness properties of the solution. The role of the intrinsic growth rate, carrying capacity, inde- pendent consumption and harvesting quota (the two last ones being closely related to the “impulsive” characteristic of the model at sampling instants) are explained and their relevance to the properties of the model are emphasized in the various given results. The concepts of best fitness, or dominance, of some species on the remaining ones, as well as the species consensus, are given and some related results are addressed. Section 3 is concerned with the Lyapunov stability properties of the equation, with special emphasis given to the case when its parameterization is time-varying. Some specific extinction conditions are investigated in Section 4, while some numerical investigation is discussed in Section 5. Finally, our conclusions end the paper. 2. The Competition Equation and Its Main Properties 2.1. The Single-Species Impulsive Beverton–Holt Equation The subsequent impulsive discrete Beverton–Holt equation for a single species, [26,27,33] extends the standard one, which has been described in the introduction: µ(k)K (k) x (k) x − ( k + 1) = , K ( k ) + ( µ ( k ) − 1) x ( k ) x ( k + 1 ) = x + ( k + 1 ) = b ( k + 1 ) x − ( k + 1 ) − d ( k + 1 ) ; ∀ k ∈ N0 (1) Appl. Sci. 2021, 11, 9020 4 of 38 with given finite initial conditions x (0) ≥ 0, where µ(k) ∈ R0+ := R+ ∪ {0} is the so-called intrinsic growth rate of the population, determined by life cycle and demographic properties (such as, for instance, species growth rate, survivorship rate, etc.) and K (k) ∈ R+ is the so-called carrying capacity, a characteristic of the habitat dependent on resource availability, temperature, humidity, etc. Typically, the constraints µ(k) > 1 and K (k) ∈ R+ for all k ∈ N0 are assumed for the Beverton–Holt equation, since if, otherwise, the carrying capacity is zero, then the habitat is incompatible with the species evolution within it; if µ(k) ≡ 1, then the species stock remains constant through time; if 0 < µ(k ) < 1, then x − (k + 1)/x (k) = µ(k )/[1 − (1 − µ(k)) x (k )/K (k )], so that two special atypical situations can happen, namely, if 1 > µ(k) > 1 − (1 − µ(k)) x (k)/K (k) for the initial samples and there is a sufficiently small population stock, which keeps the right-hand-side positivity, then the species stock grows over time, but it can tend to negative values, making the model to be not well-posed for a sufficiently extended time. If, on the contrary, µ(k) < 1 − (1 − µ(k )) x (k)/K (k) for all samples, then the population asymptotically extinguishes. On the other hand, the impulsive parameters, which can modify the stocks of pop- ulation, to the right of the sampling instants related to their left-values are the harvest- ing quota (1 − b(k)) ∈ R+ ; ∀k ∈ N0 , which depends on the authorities policy depend- ing on the spawning stock and foreseen recruitments and the independent consumption d(k )(≤ b(k ) x (k)) ∈ R; ∀k ∈ N0 that implies recruitment variations depending on un- foreseen disturbances, artificial re-population actions and local migrations “towards” or “from” the studied habitat. For easy presentation of some of the technical results, the complementary harvesting quota b(k) ∈ R0+ ; ∀k ∈ N0 is used. Some particular results in this paper are given for the case of intrinsic growth rates under unity and for negative harvesting quotas, or, equivalently, for complementary harvesting quotas exceeding unity. The physical interpretation of this situation is that a re-population with a number of adults exceeding to population stock is performed with a ratio larger than unity. Remark 1. Note that the public regulations on fishing or hunting taking into account the stocks of population imply that b(k) ≤ 1; ∀k ∈ N0 . However, if a repopulation with adults is intended, which is sometimes the case, then it can be decided to increase the above coefficient exceeding unity in the case that the population stock is low. Therefore, we admit, at the level of a general formulation, that the complementary harvesting quota is non-negative and bounded. It turns out that the independent consumption is positive, if there is a net flux of individuals from the studied habitat to outside of it, so that the population stock decreases, and it is negative, if there is a net flux entering the habitat, so that the population stock increases. Finally, it is zero if there is no interchange of population of the habitat with the environment. The independent consumption can also take into account repopulation actions that are not linked to a proportion on the existing population stock. Note that, if x (0) = 0, or if K (0) = 0 and {d(k)}∞ ∞ k =0 ≡ 0, then { x ( k )}k =0 ≡ 0. It can be argued that, although it is a common hypothesis to assume that {K (k)} ⊂ R+ and {µ(k) − 1} ⊂ R+ , it is possible, under certain extra constraints on the state of sufficient smallness, related to information about the equation parameters, that µ(k) could take values under unity. To clarify this concern, let us assume Equation (1) without impulsive effects, that is, {b(k )}∞ ∞ − k =0 ≡ 1 and { d ( k )}k =0 ≡ 0, so that x ( k ) = x ( k ); ∀ k ∈ N0 . Thus, K (k) for any given k ∈ N0 , if µ(k ) ∈ [0, 1], then x (k + 1) ≥ 0 if 0 ≤ x (k) ≤ 1− µ ( k ) , and x (k + 1) Furthermore, if µ(k) = 0, then { x ( j)}∞ K (k) is, furthermore, finite if x (k ) < 1− µ ( k ) . k ≡ 0, so that extinction happens in finite time. If {µ( j)}∞ ∞ k ≡ 1, then { x ( j )}k ≡ x ( k ), so that the population remains constant and fixed to its value at the k − th sample. Appl. Sci. 2021, 11, 9020 5 of 38 K (k) Now, let us assume that µ(k) ∈ (0, 1) and 0 < x (k) < 1− µ ( k ) , implying that x (k + 1) > K ( k +1) 0 and finite. Let us also assume that µ(k + 1) ∈ (0, 1). Then, 0 < x (k + 1) < 1−µ(k+1) for x (k + 2) > 0 and finite. One obtains µ(k)K (k) x (k) K ( k + 1) x ( k + 1) = < K (k) − (1 − µ(k)) x (k) 1 − µ ( k + 1) which implies that [µ(k)K (k)(1 − µ(k + 1)) + K (k + 1)(1 − µ(k))] x (k) < K (k)K (k + 1) K (k) so that the former non-negativity constraint x (k) < 1− µ ( k ) is amended by the stronger one K ( k ) K ( k + 1) K (k) x (k) < < µ(k )K (k)(1 − µ(k + 1)) + K (k + 1)(1 − µ(k )) 1 − µ(k) And, proceeding recursively in the same way, it follows that x (k) has a strictly decreas- ing infinite sequence of bounded non-negative upper-bounds. However, note that, if µ(k + 1) = 1, then the upper-bound remains unaltered. Thus, we have the following results. Lemma 1. Let us assume that {d(k)}∞ ∞ k =0 ≡ 0 and { b ( k )}k =0 ≡ 1 (namely, there is no impulsive effect). Let us also assume that, for some finite k ∈ N0 , {µ( j)}∞ ∞ k0 ⊂ (0, 1), { x ( j )}k0 ⊂ R0+ is bounded and either {K (k)}∞ ∞ k =k0 is non-decreasing or { K ( k )} j=k0 ⊂ (0 , 2). Then, ∃ k 1 ≥ 1 such that x (k) = 0 for k ≥ k0 + k1 (that is, there is extinction in finite time). Proof. Note that, since {µ( j)}∞ ∞ k0 ⊂ (0, 1), { x ( k )}k=0 ⊂ [0, ∞ ) is strictly decreasing if and only if   µ(k)K (k) x ( k + 1) − x ( k ) = − 1 x (k) < 0; ∀k(≥ k0 ) ∈ N0 K (k) − (1 − µ(k )) x (k) If x (k) = 0 for some finite k, then the result is already proved. Thus,  let us assume  that K (k) x (k ) > 0; ∀k(≥ k0 ) ∈ N0 . Note that x (k) > 0, if and only if x (k) < min K (k), 1−µ(k) = µ(k)K (k) K (k ), since K (k)−(1−µ(k)) x(k) < 1, that is, since (1 − µ(k)) x (k) < (1 − µ(k))K (k ), leading to x (k ) < K (k); ∀k(≥ k0 ) ∈ N0 , in order for it to be jointly strictly decreasing with the denominator of the right-hand-side being positive for all ∀k(≥ k o ) ∈ N0 . Then, either ( x (k) − K (k)) → 0 as k → ∞ , or x (k) → 0 as k → ∞ ; furthermore, x (k) < min(K (k), x (k − 1)) < . . . < inf K ( j); ∀k(≥ k0 ) ∈ N0 k ≥ j ≥0 Since {K (k)}∞ j=k0 is non-decreasing, then inf K (k) = inf K (k) ≥ inf K (k ), so k 0 ≥ k ≥0 ∞ ≥ k ≥0 ∞≥k≥k0 that   x (k) 0. Let us also ∞ assume that, for some finite k ∈ N0 , {µ( j)}k ⊂ (0, 1], card{ j(≥ k) ∈ N0 :µ( j) < 1} < ∞ and { x ( j)}∞ k ⊂ R0+ is bounded. Then, x ( k ) > 0, that is, there is no extinction in finite time and x (k ) = x (> 0); ∀k ≥ k0 (∈ N0 ) with finite k0 , so that a non-extinction equilibrium point is reached in finite time. In Lemma 2, note that the condition card{ j(≥ k ) ∈ N0 :µ( j) < 1} < ∞ implies that card{ j(≥ k) ∈ N0 :µ( j) = 1} = χ0 , that is an infinity denumerable cardinal. In other words, there is only a finite number of samples in the infinite subset of the natural numbers {k, k + 1, . . . , +∞} for which the intrinsic growth rate is strictly less than unity. This condition guarantees that, contrarily to Lemma 1, there is no extinction at the k − th sample. Additionally, it turns out that, since the intrinsic growth rate sequence is less than one for a finite number of samples, there exists a finite non-negative integer number k0 such that µ(k ) = 1: ∀k ≥ k0 implying that x (k0 + j) = x (k0 ); ∀ j ∈ N0 leading to the existence of an equilibrium point. It is addressed, in Section 5, how it is possible to achieve asymptotic extinction in the single species and in the competition impulsive Beverton–Holt equations under certain trade-offs between the independent consumption sequence and the harvesting quota, even if the intrinsic growth rate is under unity. 2.2. The Competition Impulsive Beverton–Holt Equation The class of impulsive discrete Beverton–Holt Equation of n competing species pro- posed for further analysis is a generalization of the discrete Beverton–Holt Equation (1) for the particular case of a single species, i.e., for n = 1, with a constant parameterization of in- trinsic growth rates and carrying capacities through time being expressed for n competing species as µ i Ki x i ( k ) xi− (k + 1) = n (2) Ki + ∑ j=1 c ji (µi − 1) x j (k) xi (k + 1) = bi (k + 1) xi− (k + 1) − di (k + 1); ∀i ∈ n = {1, 2, . . . , n}, ∀k ∈ N0 (3) with given finite initial conditions xi (0) ≥ 0. For some coupling coefficients cij ≥ 0; ∀i, j ∈ n with cii = 1; ∀i ∈ n. The subscripts stand for the populations and model parameters of the various competing species, while each current sample is denoted in parenthesis. The superscript “-” stands for the population values to the left of a sampling instant, while the absence of superscript indicates its value to the right of sampling instant. Note that it was assumed, for the sake of exposition simplicity, that the intrinsic growth rate, carrying capacity and inter-species coupling dynamics coefficients are constant through time for each one of the species. The impulsive control parameters are the members of the real sequences {bi (k)}∞ ∞ k =1 and { di ( k )}k =1 ; ∀i ∈ n. Those sequences may modify the levels of populations to the right of the sampling instants. If those sequences are identically zero for any population at a sampling instant, then its value after sampling remains unaltered. The combination of (2) and (3) yields   [bi (k + 1)µi Ki − di (k + 1)(µi − 1)] xi (k) − di (k + 1) Ki (k) + ∑nj(6=i)=1 c ji (k)(µi − 1) x j (k) x i ( k + 1) = ; ∀k ∈ N0 , ∀i ∈ n. (4) Ki () + ∑nj=1 c ji (µi − 1) x j (k) It is of interest to summarize, first, some basic properties of the impulsive time- invariant single Beverton–Holt equation for only one species to then interpret some parallel properties associated with a tandem of competing species. It is apparent that the equi- librium point exists also if the left and right values at the sampling time instants are non-negative and finite, even if they are distinct. In the case that they are distinct, there is an oscillatory equilibrium solution switching among those left and right values at the sampling time instants. In the case that the four parameters are potentially time-varying, Appl. Sci. 2021, 11, 9020 7 of 38 the mentioned impulsive Beverton–Holt equation for a single species (that is, n = 1) is the following one: µ(k)K (k) x (k) x − ( k + 1) = ; K ( k ) + ( µ ( k ) − 1) x ( k ) x ( k + 1) = b ( k + 1) x − ( k + 1) − d ( k + 1) that, combined, yield x ( k + 1) = b ( k + 1) x − ( k + 1) − d ( k + 1) [µ(k)b(k+1)K (k)−d(k+1)(µ(k)−1)] x (k)−d(k+1)K (k) = K (k )+(µ(k )−1) x (k ) 2.3. Equilibrium Points The subsequent result holds for the case of one single species. Theorem 1. Let us assume that either the intrinsic growth rates, carrying capacities, harvesting quotas and independent consumptions are constant parameters, or that the related parameterizing sequences converge asymptotically to constant parameters. Then, the following properties hold: (i) There exist equilibrium points (or, more properly, left and right limits of equilibrium points) q µbK − d(µ − 1) − K ± (µbK − d(µ − 1) − K )2 − 4(µ − 1)dK x1,2 = 2( µ − 1) − x1,2 + d x1,2 ≥ b with physical sense, or feasible, (i.e., being real non-negative) of the single-species impulsive Beverton–Holt equation if and only if " p 2 !# " p 2 ! µb − 1 K (µb − 1)K µb + 1 K d ∈ − x1 , min , ∪ ,∞ µ−1 µ−1 µ−1 If µ ≥ max 1, b−1 , then the above constraint becomes  " p 2 # " p 2 ! µb − 1 K µb + 1 K d ∈ − x1 , ∪ , ∞ µ−1 µ−1 (ii) If d = 0 and b = 1, then there are no impulses at the equilibrium steady-state, so that the equilibrium solution is not oscillatory between its left and right values at the sampling instants, that is, x1 = x1− = 0 (extinction) and x2 = x2− = K. (µb−1)K If d = 0, then the equilibrium points become x1 = 0 (extinction) and x2 = µ−1 > 0, provided that either µ > max 1, b−1 or µ < min 1, b−1 . If, in addition, µ > 1 and   b = µ−1 , then x1 = x2 = 0. r (µb−1)K d (iii) If d = µ −1 < 0, then there is only a non-negative equilibrium point x1 = µ −1 K>0 with either µ > 1 and d < 0 or µ < 1 and d < 0. √ µbK 1+4(1−µb)−1 (iv) If d = then there is only a non-negative equilibrium point x1 = µ −1 , 2( µ −1) K>0 − 1   with µ > 1 and b ∈ 0, µ . Let us assume that negative harvesting quotas are allowed, √ 1+ 4µb−3 then x1 = 2(1−µ) > 0, provided that b > max µ−1 , 34 µ .  Appl. Sci. 2021, 11, 9020 8 of 38 Proof. The equilibrium point to the left and right of the sampling instants satisfies x − = x +d b subject to the constraint [µb K − d(µ − 1)] x − dK x= K + ( µ − 1) x The above constraint may be rewritten as (µ − 1) x2 − [µbK − d(µ − 1) − K ] x + dk = 0 whose zeros x1 , x2 ≥ x1 are q µbK − d(µ − 1) − K ± (µbK − d(µ − 1) − K )2 − 4(µ − 1)dK x1,2 = , 2( µ − 1) which have to be real and non-negative for the equilibrium point to have a physical sense, so that two joint necessary conditions are (µb − 1)K d ≤ d0 = ; (µbK − d(µ − 1) − K )2 ≥ 4(µ − 1)dK µ−1 The second constraint can be re-arranged in the form (µ − 1)2 d2 − 2[(µb − 1) + 2]K (µ − 1)d + (µb − 1)2 K2 ≥ 0 Defining the auxiliary variable f = (µ − 1)d, the above inequality takes the form g( f ) = f 2 − 2(µb + 1)K f + (µb − 1)2 K2 ≥ 0 One obtains that the zeros of the left-hand-side of the above constraint are  q   2 2 p  f 1,2 = µb + 1 ± (µb + 1) − (µb − 1) K = µb + 1 ± 2 µb K and the constraint g( f ) ≥ 0 holds, since g( f ) = 0 is a convex parabola of zeros f 1 , f 2 ≥ f 1 , if any of the two constraints below holds p 2 p 2 µb + 1 K µb − 1 K d ≥ d2 = or d ≤ d1 = µ−1 µ−1 − x +d In order that, in addition, x1,2 ≥ 1,2b ≥ 0 (in order for the equilibrium points to the left of the sampling instants to be non-negative as well), d ≥ − x1 since 0 ≤ x1 ≤ x2 . 2 Furthermore, d ∈ [− x1 , min(d0 , d1 )] ∪ [d2 + ∞). If µ ≥ max 1, b−1 , then  p µb − 1 ≤ µb − 1. Property (i) follows directly from the above constraints. Property (ii) follows directly from the corresponding particular cases of Property (i). (µb−1)K To prove Property (iii), note that, if d = µ−1 < 0, then the equilibrium point r d satisfies the constraint (µ − 1) x2 − Kd = 0, which yields x = µ−1 K > 0 for µ > 1 and d > 0, or for µ < 1 and d < 0. The equilibrium point is unfeasible for alternative combinations of values of µ and d. µbK To prove Property (iv), note that, if d = µ−1 , then the equilibrium point constraint [µb K −d(µ−1)] x −dK |d|K (µb−1)K x= K +(µ−1) x becomes x = K +(µ−1) x for d = µ −1 < 0. Thus, the equivalent constraint (µ − 1) x2 + Kx − |d|K = 0 has, as zeros, p p −K ± K2 + 4(µ − 1)|d|K −1 ± 1 + 4(1 − µb) x1,2 = = K 2( µ − 1) 2( µ − 1) Appl. Sci. 2021, 11, 9020 9 of 38 Let us consider two cases for d < 0, namely, Case a: µb > 1 and µ < 1 so that, if b > µ−1 > 1 is not admissible, √ so that x1,2 < 0 (i.e., 1+ 4µb−3 both unfeasible since negative). If b > 1 is admissible, then x1 = 2(1−µ) > 0, provided √ − 1 3  1− 4µb−3 that b > max µ , 4 µ and x2 = 2(1−µ) < 0 (i.e., unfeasible), since µb > 1. √ −1− 1+4(1−µb) Case b: µb < 1 and µ > 1, implying b ∈ 0, µ−1 , leading to x2 =   2( µ −1) K< √ −1+ 1+4(1−µb) 0 (unfeasible) and x1 = 2( µ −1) K > 0. Therefore, only one equilibrium point is non-negative. Property (iii) has been proved.  Note that Theorem 1(iii) establishes the existence of a non-negative equilibrium point with intrinsic growth rates exceeding unity and positive harvesting quotas subject to a prescribed upper-bound. The physical interpretation is that there is a certain migration flux (either spontaneous or due to a monitored re-population intervention) of adult individuals from the habitat environment entering the habitat under study. Note that Theorem 1(iii) was also extended by first considering that Case a is non-negative with b > µ−1 > 1, with µ < 1 and d < 0. This implies that there is a negative migration flux (i.e., a supply of individuals from the environment to the habitat, either spontaneous or via intervention) and a negative harvesting quota, namely, a monitored supply of larvae to the habitat with associated complementary harvesting quota in a proportion √ greater than unity related 1− 4µb−3 to its natural stock within the habitat. In this case, x2 = 2(1− µ ) < 0 (unfeasible) and √ 1+ 4µb−3 x1 = 2(1−µ) > 0, provided that b > max µ−1 , 34 µ .  2.4. Non-Negativity and Boundedness of the Solutions The subsequent assumptions are eventually made in the sequel in some of the technical results. Assumption 1. The carrying capacities, intrinsic growth rates and inter-species coupling dynamics are constant, i.e., Ki > 0, µi > 1, cii = 1, cij ≥ 0; ∀i, j(6= i ) ∈ n. Assumption 2. The complementary harvesting quotas and independent consumption sequences satisfy the upper-bounding constraints 0 ≤ bi (k ) ≤ bi < +∞; −∞ < di (k) ≤ bi (k) xi− (k ); ∀i ∈ n, ∀k ∈ N0 . It can be pointed out that Assumption 1 can be extended directly to the case that the constant parameters are sequences depending on time. It can be also generalized for µi ≤ 1 Ki n under the additional constraint 1− µi ≥ ∑ j=1 c ji x j ( k ); ∀i, j (6 = i ) ∈ n or for its extension to the case when carrying capacities, intrinsic growth rates and inter-species coupling dynamics are sequences depending on time. The interpretation is that, if the intrinsic growth rate is under unity, then the coupling sequences from the remaining populations to the i − th have to be small enough related to the ratio Ki (k)/(1 − µi (k )). Note that the populations stocks to the left and right limits of the sampling instants have to be non-negative for coherency of the model evolution. The constraints of Assumption 2 guarantee that the populations to the left limits of the sampling instants are non-negative if those to the right limits are non-negative. The following result on the non-negativity of the solution is directly obtained from (2)–(3) and Assumptions 1–2. Proposition 1. Under Assumptions 1–2, for all i ∈ n and any given k ∈ N0 , if xi (k) ≥ 0, then xi− (k + 1) ≥ 0 and xi (k + 1) ≥ 0. If xi (k) > 0, then xi− (k + 1) > 0 and, also, xi (k + 1) > 0 if, furthermore, di (k + 1) < bi (k + 1) xi− (k + 1). Appl. Sci. 2021, 11, 9020 10 of 38 Assumptions 1–2 also lead to the subsequent result on the dominance (or best fitness) of a certain population on another one for all samples. Proposition 2. Let us assume that, for any given i, j(6= i ) ∈ n, µi ≥ µ j , Ki ≥ K j = ∞ ∞   µ i ( µ i −1) max cc`i   K j , bi (k) − b j (k ) k=1 ⊂ R0+ and d j (k ) − di (k) k=1 ⊂ R0+ and that `∈n ` j µ j ( µ j −1 ) n o ∞ Assumptions 1–2 hold. Thus, if xi (0) ≥ x j (0), then xi− (k) − x − j (k) ⊂ R0+ and k =0  ∞ x i ( k ) − x j ( k ) k = 0 ⊂ R0 + . The above conditions may be relaxed as follows: n o∞ µ If Ki ≥ K j , µi ≥ µ j and xi (0) ≥ µj x j (0), then xi− (k ) − x − j ( k ) ⊂ R0+ and i k =1 n o∞ xi− (k) − x − j (k) ⊂ R0+ , provided that d j (k) − di (k) ≥ b j (k) x − − j ( k ) − bi ( k ) xi ( k ); ∀ k ∈ N. k =1 Proof. Note that µ i x i (0) xi− (1) = 1+∑n`=1 c`i (µi −1)Ki−1 x` (0) µ j x j (0) (5) ≥ x− j (1) = 1+∑n`=1 c` j ( µ j −1 ) K− 1 j x ` (0) if xi (0) ≥ x j (0), µi ≥ µ j and   1 n 0≤ ∑ `=1 c`i (µi − 1)Ki−1 − c` j µ j − 1 K − j x ` (0) ≤ 0 which holds if Ki ≥ K j and µi xi (0) ≥ µ j x j (0), which leads to xi− (1) ≥ x −j (1). Furthermore, xi (1) ≥ x j (1) since xi (0) ≥ x j (0), bi (1) − b j (1) ≥ 0 and d j (1) − di (1) ≥ 0. Then, one o∞ ∞ n obtains recursively that xi− (k) − x −  j ( k ) ⊂ R0+ and xi (k) − x j (k) k=0 ⊂ R0+ , as k =0 claimed in the first part of the proposition. The second part is direct, since, under the µ assumptions that Ki ≥ K j and xi (0) ≥ µj x j (0), one obtains µi Ki xi (0) ≥ µ j K j x j (0) (not i necessarily under the joint constraints xi (0) ≥ x j (0), µi ≥ µ j and Ki ≥ K j ). This implies that xi− (1) ≥ x − − − j (1), which, together with d j (1) − di (1) ≥ b j (1) x j (1) − bi (1) xi (1) (not necessarily with d j (1) ≥ d1 (1)), implies that xi (1) ≥ x j (1). The result for k ≥ 1 follows recursively.  Remark 2. It follows, from Proposition 2, that a population i ∈ n is dominant on the remaining ones, or, in other words, it has the best fitness, for all j(6= i ) ∈ n and the respective  left and right  sam- 1 −  pled values if Ki ≥ max K j , xi (0) ≥ µ max µ j x j (0) and min d j (k) − b j (k) x j (k) ≥ j(6=i )∈n i j(6=i )∈n j(6=i )∈n di (k) − bi (k ) xi− (k); ∀k ∈ N. It turns out, furthermore, that an increase in an individual carrying capacity facilitates the increase in the corresponding species numbers, if all the remaining parameters remaining unchanged. The above constraints generalize Proposition 2, although at the expense of involving conditions on the left amounts of the population stocks. It can be pointed out that, eventually, there could exist more than one species with an identical best fitness within the set of competing species. The following result is now proved on non-negativity, boundedness and eventual extinction of the species populations. Theorem 2. The following properties hold: (i) If ∞ > xi (0) > 0; ∀i ∈ n, Assumption 1 holds, {bi (k)}0∞ ≡ 1 and {di (k)}0∞ ≡ 0; ∀i ∈ n; then, { xi (k )}∞ ∞ k =0 is non-negative and bounded; ∀i ∈ n. Furthermore, if { xi ( k )}0 → 0 (i.e., asymptotic extinction of the ith species), then lim inf min Kie (k) = 0. k→∞ 0≤ k < ∞ Appl. Sci. 2021, 11, 9020 11 of 38 (ii) If ∞ > xi (0) > 0; ∀i ∈ n and Assumptions 1 and 2 hold with {di (k)}∞ k =0 ⊂ R0+ , then and it is a bounded non-negative sequence if ∏∞ k =0 [ µ i b i ( k )] = 0; ∀ i ∈ n. This property holds, in particular, if bi (k ) > µi−1 ; ∀k ∈ N0 . If, furthermore, 0 ≤ di (k ) < bi (k) xi− (k); ∀k ∈ N0 , then extinction in finite time does not occur for any i ∈ n, so that extinction only happens asymptotically and { xi (k)}∞ ∞ k =0 is a bounded non-negative sequence if ∏k =0 [ µi bi ( k )] ≤ A < j + a −1 ∞. In particular, if ∏k= j [bi (k)] ≤ µi− a ; ∀k (≥ j0 ) ∈ Z0+ and some finite j0 ∈ N0 and a ∈ N, then µ i Ki x i ( k ) xi− (k + 1) = n (6) Ki + ∑ j=1 c ji (µi − 1) x j (k) If Assumptions 1 and 2 hold with {di (k)}∞ ∞ k=0 ⊂ R, then { xi ( k )}k =0 → 0 if ∞ ∏k=0 [µi bi (k) − di (k + 1)] = 0. (iii) If Assumptions 1 and 2 hold, then { xi (k)}∞ k =0 is bounded; ∀i ∈ n. Proof. One obtains that (2) is equivalent to µi xi ( k ) µi xi ( k ) xi− (k + 1) = µ i −1 n = µ i −1 ; ∀ k ∈ N0 (7) 1+ Ki ∑ j=1 c ji x j ( k ) 1+ x (k) Kie (k ) i where the equivalent i − th carrying capacity Kie is given by Ki x i ( k ) Ki 0 ≤ Kie (k) = n = ≤ Ki ; ∀k ∈ N0 (8) ∑ j=1 c ji x j (k) 1 + ∑nj(6=i)=1 c ji xj (k) x (k) i since cii = 1 and cij ≥ 0; ∀ j(6= i ) ∈ n. By combining (7) and (8), one obtains µi xi ( k ) µi xi ( k ) µ i Ki x i ( k ) xi− (k + 1) = µ −1 ≤ µ i −1 = Ki +(µi −1) xi (k ) ; 1+ K i ( k ) x i ( k ) 1+ Ki x i ( k ) (9) ie ∀i ∈ n, ∀k ∈ N0 If Assumption 1 holds, {bi (k)}0∞ ≡ 1 and {di (k)}0∞ ≡ 0, then xi− (k) = xi (k); ∀i ∈ n, ∀k ∈ N0 . Let us define yi (k) = 1/ xi− (k) ; ∀i ∈ n, ∀k ∈ N0 . One obtains, from recursive calculations in (8) using (7), j−k µ i −1 µ i −1 µ i −1 k µi 0< µ i Ki ≤ µi max Kie (k ) ≤ µi ∑ j=0 Kie ( j) 0≤ k < ∞ µ i −1 ≤ 1 x i ( k +1) = yi (k + 1) = yi− (k + 1) = µi−1 yi (k) + µi Kie (k ) j−k −(k+1) µ i −1 k µi = µi y i (0) + µi ∑ j=0 Kie ( j) µ i −1 k j−k µi min Kie ( j) ∑ j=0 i ≤ y i (0) + µ 0≤ j ≤ k j−k (10) 1 µ i −1 ∞ µi min Kie ( j) ∑ j=0 ≤ x i (0) + µi 0≤ j ≤ k −(k+1) 1 µ i −1 1− µ i ≤ x i (0) + µi min Kie ( j) 1−µ−1 0≤ j ≤ ∞ i 1 1 ≤ x i (0) + min Kie (k) ; 0≤ k < ∞ ∀i ∈ n, ∀k ∈ N0 since Kie (k) ≤ Ki < ∞. Since ∞ > xi (0) > 0 and > 0; ∀i ∈ n, ∀k ∈ N0 , then 1 x i ( k +1) { xi (k)}0∞ is bounded; ∀i ∈ n. Note that { xi (k)}0∞ is a non-negative sequence for any given i ∈ n by construction and, furthermore, a positive sequence, unless lim inf min Kie (k) = 0, k→∞ 0≤ k < ∞ which would eventually imply, from (10), that {1/xi (k )}0∞ → ∞ , so that { xi (k)}0∞ → 0 . ∞ Since it has been already proved (Property (i)) that x j (k) 0 is bounded for any j ∈ n, Appl. Sci. 2021, 11, 9020 12 of 38 ∞ → ∞ , which would translate into lim inf min Kie (k) = 0 from (8),  then x j (k)/xi (k) 0 k→∞ 0≤ k < ∞ if and only if { xi (k)}0∞ → 0 . Property (i) has been proved. From Assumptions 1–2, with x ( k +1) non-negative independent consumption, one directly obtains that ix (k) ≤ bi (k + 1)µi , i which leads to the first part of Property (ii). If Assumptions 1 and 2 hold with non- necessarily non-negative independent consumption sequence, then { xi (k)}∞ k=0 → 0 if ∏∞ [ µ k =0 i ib ( k ) − d i ( k + 1 )] = 0, which proves the second part of Property (ii). Property (iii) is proved by contradiction arguments. Let us assume that Assumptions 1 and 2 hold and { xi (k )}∞ k =0 is unbounded. Thus, there is a subsequence { xi ( nk )}k =0 ⊂ ∞ ∞ { xi (k)}, with nk+1 − nk < +∞, which is strictly increasing, so that { xi (nk )}k=0 → +∞ and then, from (7), x − ( k + 1) µi Kie (k ) 0≤ i ≤ ; ∀ k ∈ N0 xi ( k ) ( µ i − 1) x i ( k ) xi− (nk +1) then lim = 0. Furthermore, k → ∞ xi ( n k )    xi− (nk+1 ) xi− (nk+1 ) x i ( n k +1 −1 ) xi ( n k ) = x i ( n k +1 −1 ) xi− (nk+1 −1)      xi− (nk+1 −1) x i ( n k +1 −2 ) xi− (nk +1) (11) × x i ( n k +1 −2 ) xi− (nk+1 −2) ... xi ( n k ) x − ( n +1) ≤ C (nk , j) i x (nk ) i k where C (nk , j) is a positive real constant which is finite, since it is bounded to a finite number of iterations between nk and nk+1 , since nk+1 − nk < +∞ and since one has, from Assumption 2, that xi (k + 1) = bi (k + 1) xi− (k + 1) − di (k + 1) ≤ Mi < +∞; ∀i ∈ n; ∀k ∈ N0 and x i ( k + 1) d ( k + 1) M 0≤ − = bi (k + 1) − −i ≤ − i < +∞ (12) x i ( k + 1) x i ( k + 1) x i ( k + 1) x i ( k +1) if xi− (k + 1) 6= 0 and = +1 < +∞ if xi− (k + 1) = 0, still from Assumption 2. xi− (k +1) = 0 0  − ∞ x (n ) x i ( n k +1) Then, from (11) and (12), one obtains lim x (kn+1) = 0, since x (n ) is bounded and k→∞ i k i k k =0 { xi (nk )}∞ k =0 is strictly increasing from the contradictionargument  hypothesis. Thus, there   ∞ x i m n k +1 is a subsequence xi mnk k=0 ⊂ { xi (nk )} such that < 1, so that one concludes, xi ( m nk )   ∞ from D´Alembert criterion for convergence of sequences, that xi mnk k=0 → 0 . As a result, { xi (k)}∞k =0 is bounded. Property (iii) has been proved.  Note that Theorem 2 (i) guarantees the boundedness of all the populations under Assumption 1, provided that the harvesting quota and the independent consumption are zero (so that there are no impulsive effects at sampling instants), while Theorem 2(iii) guarantees such a boundedness under Assumptions 1 and 2. In the following, the limits of the convergent sequences are denoted with the same notation as the sequences of elements by suppressing the sample argument “k” for the sake of notation simplicity. Theorem 3. Let us assume that all the populations converge to equilibrium points, that is, { xi (k)}0∞ → xi . Then, the following properties hold: (i) Let us assume that Assumptions 1 and 2 hold and that {bi (k)}0∞ → bi > µi−1 ,   ∞ {di (k)}0∞ → di ∈ 0, di , c ji (k) 0 → c ji ∈ 0, c ji for some di ∈ R0+ , c ji ∈ R0+ ; h i    Appl. Sci. 2021, 11, 9020 13 of 38 ∀ j(6= i ), i ∈ n. If di and c ji ; ∀ j(6= i ), i ∈ n are sufficiently small compared to Ki (bi µi − 1), then there are two equilibrium points xi ≥ 0; ∀i ∈ n. (ii) Let us assume that the quadruple (Ki , µi , bi , di ) satisfies the constraints Bi = [Ki (bi µi − 1) − (µi − 1)di ] ≥ 0 Ci = (Ki (bi µi − 1) − (µi − 1)di )2 − 4(µi − 1)di Ki ≥ 0 Then, one of each i − th equilibrium point is feasible, that is, xi ≥ 0 under the subsequent stipulations:  √  B + B2 − C (1) Either the coupling constraint ∑nj(6=i)=1 c ji x j ∈ i µ −1i i , ∞ from the remaining i equilibrium points and, furthermore, any of the three conditions below: K i ( bi µ i − 1 ) C1) bi ≥ µi−1 and 0 ≤ di ≤ µi − 1 Ki (bi µi − 1) 4Ki  i   C2) µi > 1, bi ∈ µi−1 , , di ∈ 5µi−1 , µi − 1 µi − 1 K i ( bi µ i − 1 )    −1 −1 i 4Ki C3) µi < 1, bi ∈ µi , 5µi , di ∈ − , − ; 1 − µi 1 − µi (2) Or, the coupling constraint ∑nj(6=i)=1 c ji x j ∈ [0, ∞) from the remaining equilibrium points together with sgn(di ) = sgn(µi − 1) (which excludes the values µi = 1 and di = 0). Proof. Take any i ∈ n. One obtains, from (2)–(3), by taking limits, that µ i Ki x i xi− = (13) Ki + ∑nj=1 c ji (µi − 1) x j bi µ i K i x i xi = bi xi− − di = Ki +∑nj=1 c ji (µi −1) x j − di   (14) (bi µi Ki −di (µi −1)) xi −di Ki +(µi −1)∑nj(6=i)=1 c ji x j = Ki +∑nj=1 c ji (µi −1) x j then   di Ki + (µi − 1)∑nj(6=i)=1 c ji x j ! (bi µi Ki − di (µi − 1)) − 1 xi = (15) Ki + ∑nj=1 c ji (µi − 1) x j Ki + ∑nj=1 c ji (µi − 1) x j so that   di Ki +(µi −1)∑nj(6=i)=1 c ji x j xi = bi µi Ki −di (µi −1)−Ki −∑nj=1 c ji (µi −1) x j   (16) di Ki +(µi −1)∑nj(6=i)=1 c ji x j = h i (µi −1) (Ki −di − xi )−∑nj(6=i)=1 c ji x j +(bi −1)Ki µi by noting that the denominator of the above second identity balances that of the first one. Re-arranging terms in the above Equation (16) leads to (µi − 1) xi2 − b̂i xi + ĉi = 0 (17) where b̂i = b̂i ( x1 , . . . , xi−1 , xi+1 , ·· · , xn ) = (µi − 1) Ki − di − ∑nj(6=i)=1 c ji x j + (bi − 1)Ki µi (18)   = Ki (bi µi − 1) − (µi − 1) di + ∑nj(6=i)=1 c ji x j Appl. Sci. 2021, 11, 9020 14 of 38   n ĉi = ĉi ( x1 , . . . , xi−1 , xi+1 , · · · , xn ) = di Ki + (µi − 1)∑ j(6=i)=1 c ji x j (19) It is obvious, from the above equations, that if bi > 1/µi and di = 0 and c ji = 0; ∀ j(6= i ), then b̂i = Ki (bi µi − 1) and ĉi = 0; ∀i ∈ n and one obtains the equilibrium points K ( b µ −1) xi = i µi −i1 > 0 and xi = 0 as zeros of (17). It turns out that, in the general case, the i zeros of (17), namely, q b̂i ± b̂i2 − 4(µi − 1)ĉi xi = (20) 2( µ i − 1) are both non-negative and real if b̂i > 0 and b̂i2 ≥ 4(µi − 1)ĉi which hold for each i ∈ n if di and c ji ; ∀ j(6= i ) are sufficiently small in view of (18) and (19), such that ĉi is sufficiently small, compared to Ki (bi µi − 1) (that is, sufficiently small related to the first additive term of b̂i ) from the property of the continuity of the zeros of a polynomial with respect to its coefficients. Property (i) has been proved. Property (ii) is now proved by guaranteeing the validity of b̂i2 ≥ 4(µi − 1)ĉi depending on the coupling equilibrium dynamics from the remaining species j(6= i ) ∈ n to the i − th one. First, let us define the auxiliary amount ∆i = (µi − 1)∑nj(6=i)=1 c ji x j so that, from (18) and (19), b̂i = (µi − 1)(Ki − di ) − ∆i + (bi − 1)Ki µi = Ki (bi µi − 1) − (µi − 1)di − ∆i ĉi = di (Ki + ∆i ) then g(∆i ) = b̂i2 − 4(µi − 1)ĉi = (Ki (bi µi − 1) − (µi − 1)di − ∆i )2 − 4(µi − 1)(di (Ki + ∆i )) = ∆2i − 2[Ki (bi µi − 1) − (µi − 1)di ]∆i + (Ki (bi µi − 1) − (µi − 1)di )2 − 4(µi − 1)di Ki − 4(µi − 1)di ∆i (21) = ∆2i − 2[Ki (bi µi − 1) − (µi − 1)di ]∆i + (Ki (bi µi − 1) − (µi − 1)di )2 − 4(µi − 1)di Ki = ∆2i − 2Bi ∆i + Ci ≥ 0 holds, where Bi = [Ki (bi µi − 1) − (µi − 1)di ] (22) Ci = (Ki (bi µi − 1) − (µi − 1)di )2 − 4(µi − 1)di Ki = Bi2 − 4(µi − 1)di Ki (23) Let us consider the subsequent two cases: Case 1 (g(∆i ) = 0 has real zeros ∆i1,i2 ). For the above q constraint (21) to be well-posed q with real zeros ∆i1,i2 of g(∆i ) = 0, with ∆i2 = Bi + Bi2 − Ci ≥ ∆i1 = Bi − Bi2 − Ci (according to the definition of ∆i ), so that g(∆i ) ≥ 0 if and only if ∑nj(6=i)=1 c ji x j ∈   √   Bi + Bi2 −Ci max 0, µ −1 , ∞ , a necessary condition being that Bi2 ≥ Ci , equivalently, i (µi − 1)di ≥ 0. Provided that such a necessary condition holds, sufficient conditions for ∆i2 = ≥ 0 are as follows: K ( b µ −1) (a) Bi ≥ 0, equivalently, 0 ≤ di ≤ i µi −i1 , which, combined with the above necessary i condition, also implies that bi ≥ µi−1 ; (b) Bi < 0 together with Ci ≤ 0. These two joint sufficiency-type conditions are K i ( bi µ i − 1 ) Bi2 respectively equivalent to di > µ i −1 and the implicit constraint di ≥ di = 4( µ i −1) Ki = (λi −µi +1)2 d2i K ( b µ −1) ,obtained after defining the auxiliary parameter λi = i id i 4( µ i −1) Ki for di 6= 0. Note i that the particular case di = 0 requires, for feasibility, that Bi ≥ 0 (from the necessary condition Bi2 ≥ Ci ), so that it is included in the sufficient condition Bi ≥ 0 (or di ≤ K i ( bi µ i − 1 ) µi −1 ), that is, it is included in (a). Since the given implicit constraint has to be checked ( λ i − µ i +1)2 d i for di 6= 0, it becomes identical to 1 ≥ 4( µ i −1) Ki , or equivalently identical to 4( µ i − 1) Ki 4( µ i − 1) Ki di ≤ di = 2 = ( λ i − µ i + 1) λ2i + ( µ i − 1)2 − 2( µ i − 1) λ i Appl. Sci. 2021, 11, 9020 15 of 38 µ i −1 Note that (µi − 1)λi = d i K i ( bi µ i − 1), so that, under the given necessary conditions, sgn[(µi − 1)λi ] = sign[Ki (bi µi − 1)] = sgn(bi µi − 1), which is +1 if bi > µi−1 and −1 if 4( µ i −1) Ki bi < µi−1 . Thus, if bi ≥ µi−1 , then di ≤ di if di ≤ , which is, in turn, guaranteed if λ2i +(µi −1)2  i K i ( bi µ i − 1 ) di ≤ µ4K−i1 . As a result, the sufficient condition (b) is guaranteed if di ∈ µ i −1 , 4Ki µi −1 , if i i bi ∈ µi−1 , 5µi−1 . 4( µ i −1) Ki 4Ki On the other hand, if bi < µi−1 , then di = λ2i +|µi −1|(|µi −1|+2|λi |) ≤ µ i −1 if, furthermore,  i K i ( bi µ i − 1 ) 4Ki µi > 1, implying again that the condition (b) is guaranteed if di ∈ µ −1 , µ −1 and i i 4| µ −1| K di = − λ2 +|µ −1|(|i µ −1i|+2|λ |) ≥ − 14K i −µ = 4Ki µ i −1 if µi < 1. i i i i i Note that, in both, Case 2 (g(∆i ) = 0 has complex conjugate zeros ∆i1, i2 ). This condition holds if and only if Bi2 < Ci (which also trivially requires that Ci > 0), which, in view of (23), requires that (µi − 1)di < 0, which holds if and only if sgn(di ) = sgn(µi − 1). Since g(0) = Ci > 0 and ∆i1, i2 are complex conjugate zeros, g(∆i ) > 0 in all its definition domain, that is, if ∑nj(6=i)=1 c ji x j ∈ [0, ∞). Property (ii) follows from Cases 1 and 2.  Note, from Equation (20), that both equilibrium points of the i − th species are zero if di = c ji = 0; ∀ j(6= i ) ∈ n and bi = µi−1 ; ∀i ∈ n. Excluding this trivial case, also note from Equation (20) that, provided that the radicand is non-negative, a necessary condition for at least one of the equilibrium points to be non-negative is that either b̂i > 0 or b̂i < 0 and ĉi ≥ 0. The first case implies that n K i ( bi µ i − 1 ) di + ∑ j(6=i)=1 c ji x j < µi − 1 The second case implies that µi < 1 and n K i ( bi µ i − 1 ) di + ∑ j(6=i)=1 c ji x j > ; µi − 1   n di Ki − (1 − µi )∑ j(6=i)=1 c ji x j ≥ 0 K1 while the second constraint above is fulfilled with di ≥ 0 and ∑nj(6=i)=1 c ji x j ≤ 1− µ i , or K1 di ≤ 0 and ∑nj(6=i)=1 c ji x j ≥ 1− µ i . Example 1 (Positive equilibrium points for n = 2 and n ≥ 3 with zero independent consumptions). The allocation and positivity of the equilibrium points is easy to deal with in an analytical way from a linear algebraic system in the absence of independent consumptions at the equilibrium state, as it is now discussed. If n = 2 and d1 = d2 = 0, then µ i bi K i x i xi =  for i, j(6= i ) = 1, 2 (24) Ki + (µi − 1) xi + c ji x j which leads, for xi > 0 and for i = 1, 2, denoted also as x = ( x1 , x2 ) T  0, to K i ( bi µ i − 1 ) xi + c ji x j = for i, j(6= i ) = 1, 2 (25) µi − 1 Appl. Sci. 2021, 11, 9020 16 of 38 leading to         −1 K1 (b1 µ1 −1)   K1 (b1 µ1 −1) x1 1 c21 µ1 −1 1 1 −c21 µ1 −1 =  K2 (b2 µ2 −1) = 1−c12 c21  K2 (b2 µ2 −1)  x2 c12 1 −c12 1 µ2 −1  µ2 −1  K1 (b1 µ1 −1) c21 K2 (b2 µ2 −1) (26) 1 µ1 −1 − µ2 −1 =   0 1−c12 c21 K2 (b2 µ2 −1) c12 K1 (b1 µ1 −1)  µ2 −1 − µ1 −1     provided that Ki > 0 for i = 1, 2 and either µi > max 1, bi−1 , or µi < min 1, bi−1 , for i = 1, 2 and c12 and c21 are sufficiently small according to the relations  K j ( µ i − 1) b j µ j − 1 cij <  ; i, j(6= i ) = 1, 2 (27) K i µ j − 1 ( bi µ i − 1 ) The equilibrium consensus, that is, the equilibrium points of both species are identical, x = x1 = x2 = 1−ca c > 0 for some prefixed a > 0 is achieved for the coupling coefficients 12 21 satisfying K i ( bi µ i − 1 ) µj − 1   cij = −a  ∈ [0, 1); i, j(6= i ) = 1, 2 (28) µi − 1 Kj bj µ j − 1 so that the parameter a > 0 has to satisfy, furthermore,  K i ( bi µ i − 1 ) K j bj µ j − 1 −a< ; i, j(6= i ) = 1, 2 (29) µi − 1 µj − 1 or  ! K j ( b j µ j −1 ) K i ( bi µ i − 1 ) K i ( bi µ i − 1 ) K j ( b j µ j −1 ) a > max 0, max µ j −1 − µ i −1 , µ i −1 − µ j −1 1≤i,j(6=i )≤2 (30)   K2 (b2 µ2 −1) K1 (b1 µ1 −1) = max 0, µ2 −1 − µ1 −1 The generalization to n ≥ 3 species, if di = 0; ∀i ∈ N0 , leads to a vector of equilibrium points given by   K1 (b1 µ1 −1)  −1 µ − 1 1 c21 · · ·    1 x1 cn1   ..    x2   c12 1 c32 · · · cn2   .  x=  =   .  ···  ··· ··· ··· ···    ..   xn c1n c2n · · · 1  K n ( bn µ n − 1 )  µ n −1   (31) K1 (b1 µ1 −1) µ1 −1 ..     −1  .  = ( In + o (ε))   ..     .   K n ( bn µ n − 1 ) µ n −1   provided that Ki > 0, µi > max 1, bi−1 ; ∀i ∈ n with cij ≥ 0; ∀i, j(6= i ) ∈ 3 and max cij ≤ ε for some ε ∈ R0+ so that, if ε is sufficiently small, then the inverse 1≤i,j(6=i )≤n of ( In + o (ε)) exists and x  0, i.e., all its components, that is, the equilibrium states of all the species, are positive. The following technical result, which establishes equilibrium point upper-bounds when one of the species in the competing set is dominant to the left and right samples, holds. Appl. Sci. 2021, 11, 9020 17 of 38 Proposition 3. Let us define the population inverses yi− = 1/xi− and yi = 1/xi at the equilibrium points and let us assume that the i − th species is the dominant one in the set of n competing species in the set that has the best fitness. Then, the following equilibrium constraints hold:   − ! n b − d y + d b µ − 1 µ − 1 1 Ki j(6=∑ i  i i i i yi− ≤ + i c ji i −  + −  (32) µ i bi − 1 Ki i )=1 b i b i − d i y i b i − d i y i bi −di yi− +di − i b −d y− +d i i i yi = ≤ y bi (bi −di yi− ) i(bi −di yi− ) (µi bi −1)     (33) µ i −1 µ i −1 n bi −di yi− +di 1 × + Ki ∑ j(6=i)=1 c ji + b −d y− Ki bi (bi −di yi− ) i i i If di ≤ 0 with |di | ≤ (1 − bi ) xi− , then the following further upper-bounds for (32) and (33) are obtained:     µ −1 µ −1 (1−|d |)b−1 yi− ≤ µ bbi−1 iK + iK ∑nj(6=i)=1 c ji b +|di |y−i + 1 bi +|di |yi− i i i i i i i c h  i ≤ µ bbi−1 µiK−i 1 + µiK−i 1 ∑nj(6=i)=1 b1 1 + bji i i i i   µ i −1 µ i −1 n 1 1 yi ≤ K + K ∑ j(6=i)=1 b c ji + b i i i i Proof. So that one obtains, from (2)–(3) the following equilibrium constraints, since cii = 1; ∀i ∈ n, Ki + ∑nj=1 c ji (µi −1) x j yi− = µ i Ki x i (34) µ i −1 µ i −1 n = µi−1 yi + µ i Ki + −1 µi Ki ∑ j(6=i )=1 c ji yi y j 1 yi− yi = = = bi−1 yi− + Yi (35) bi xi− − di bi − di yi− where di yi− ! 1 1 Yi = − − b yi− = (36) bi bi − di yi−  bi − d i y i i Then, from (36) into (35), one obtains bi − di yi− + di − ! di yi = bi−1 + −  yi− = y (37) bi bi − d i y i bi bi − di yi− i Which, replaced in (34), yields    di µ −1 µ −1 1 − µi−1 bi−1 + yi− = µi K + µi K ∑nj(6=i)=1 c ji yi y− 1 bi (bi −di yi− ) i i i i j   (38) = µµii−K1i + µµii−K1i ∑nj(6=i)=1 c ji y− 1 b −1 + di yi− j i bi (bi −di yi− ) and then ! !! di µ −1 n di µi − 1 µi Ki ∑ j(6=i)=1 ji j 1 − µi−1 bi−1 + − i c y−1 bi−1 + yi− = (39) bi bi − di yi− bi bi − di yi−   µ i Ki Note, from (3), that xi− ≥ max(0, di /bi ) (equivalently, yi− ≤ bi /di ) in order for the i − th equilibrium point to be non-negative to the left of the sampling instants, provided that it is non-negative to the right of such sampling instants. Therefore, the last two Appl. Sci. 2021, 11, 9020 18 of 38 additive terms on the left-hand-side of the above equation are necessarily negative under that stipulation. Let us assume that the i − th species is dominant so that xi ≥ max x j , so that 1≤ j(6=i ≤n) yi /y j ≤ 1; ∀ j(6= i ) ∈ n. Then,    di 1 − µi−1 bi−1 + yi− bi (bi −di yi− )   (40) µ i −1 µ i −1 n di ≤ + µi Ki ∑ j(6=i )=1 c ji bi−1 + µ i Ki bi (bi −di yi− ) and, since (di /bi )yi− ≤ 1,   µ i bi − 1 − µ i −1 µ i −1 n di µi−1 di µ i bi y i ≤ + µi Ki ∑ j(6=i )=1 c ji bi−1 + + y− µ i Ki bi (bi −di yi− ) bi (bi −di yi− ) i   (41) µ i −1 µ i −1 n bi −di yi− +di 1 ≤ µ i Ki + µi Ki ∑ j(6=i )=1 c ji bi (bi −di yi− ) + µi (bi −di yi− ) which is identical to (32); then, (33) is direct for the equilibrium values at the right of the sampling instants. The further two upper-bounds if di ≤ 0 inequalities follow directly from (32)–(33) and the proof is completed.  The next definition relies on the concept of equilibrium consensus among the species in the sense that all of them reach the same equilibrium points. Definition 1. It is said that the n(≥ 2) competing species set has a left/right equilibrium consensus if there exist at least two, non-necessarily distinct, non-negative real numbers x − and x such that xi− = x − and xi = x; ∀i ∈ n, that is, the equilibrium point is common to all the species. Some sufficient conditions for equilibrium consensus are found in the subsequent result. Proposition 4. The following parametric constraints hold: (i) A right equilibrium consensus for a given common population x is reached under the indepen- dent consumptions   (bi µi − 1)Ki − (µi − 1) ∑nj=1 c ji x di =   x; ∀i ∈ n (42) Ki + (µi − 1) ∑nj=1 c ji x (ii) A right equilibrium consensus for a common equilibrium population of value ( bi µ i − 1 ) K i x=   (43) (µi − 1) ∑nj=1 c ji     exists if di = 0; ∀i ∈ n with either µi > max 1, bi−1 or µi < min 1, bi−1 for each i ∈ n for x > 0, or, jointly, µi = bi−1 6= 1; ∀i ∈ n for x = 0 (that is, “extinction” consensus). (iii) Let us assume that the competing species set has an equilibrium consensus and that there are no impulses at the equilibrium points, that is, bi = 1 and di = 0; ∀i ∈ n. Then, the equivalent carrying capacities Kie = Ke ; ∀i ∈ n are identical for all the species, the carrying capacities are Ki = Ke ∑nj=1 c ji ; ∀i ∈ n and the nonzero common equilibrium point coincides with such a common equivalent carrying capacity, i.e., Ki x = x− = n = Kie = Ke ; ∀i ∈ n. ∑ j=1 c ji Appl. Sci. 2021, 11, 9020 19 of 38 Proof. One obtains, from (4), for xi = x; ∀i ∈ n that   (µi − 1) ∑nj=1 c ji x2 h  i (44) = (bi µi − 1)Ki − di (µi − 1) ∑nj=1 c ji x − di Ki ; ∀i ∈ n which proves Property (i). Property (ii) follows directly from Property (i), by making (bi µi −1)Ki  di = 0; ∀i ∈ n, which leads to x = n ; ∀i ∈ n. Since c ji ≥ 0; ∀i, j ∈ n implies (µi −1) ∑ j=1 c ji Ki µ i −1 that ∑nj=1 c ji = ≥ 0, it follows that a necessary condition of fulfillment of those bi µ i − 1 x     constraints for each i ∈ n is that either µi ≥ max 1, bi−1 or µi ≤ min 1, bi−1 , or, jointly, µi = bi−1 6= 1; ∀i ∈ n. To prove Property (iii), note that, since bi = 1 and di = 0; ∀i ∈ n, one has, at the equilibrium points, µ i Ki x x− = x = (45) Ki + ∑nj=1 c ji (µi − 1) x whose solutions are x = 0 and, after direct calculation for x 6= 0, one obtains Ki x− = x = n = Kie = Ke ; ∀i ∈ n. (46) ∑ j=1 c ji  Let us assume, for generalization purposes, that the intrinsic growth rates, the carrying capacities and inter-species coupling coefficients might also be time-varying, in general. Thus, one obtains, from the resulting generalized (2) and (3) and, provided that xi− (k + 1) 6= 0 and λ ji (k) = x j (k)/xi (k), that yi− (k + 1) = 1 xi− (k +1) (47) µi (k )−1 µi (k )−1 = µi−1 (k)yi (k) + Ki ( k ) µ i ( k ) + n Ki (k )µi (k) ∑ j(6=)=1 ji c (k)λ ji (k) yi (k + 1) = ρi−1 (k + 1)yi− (k + 1) (48) for any i, j(6= i ) ∈ n, where ρi (k + 1) = xi (k + 1)/xi− (k + 1) = bi (k + 1) − di (k + 1)yi− (k + 1) > 0 (49) for non-negativity and non-extinction in finite-time of the solution sequence, so that di (k + 1) < bi (k + 1) xi− (k + 1). Then, it follows, from (47) and (48), that yi (k + 1) = (ρi (k + 1)µi (k))−1 yi (k ) + µi (k )−1 Ki ( k ) µ i ( k ) ρ i ( k +1) (50) + K (k)µµi ((kk)− 1 n )ρ (k+1) ∑ j(6=)=1 ji ( ) ji ( ) c k λ k i i i for any i ∈ n. It can be proved that xi (k + 1) > xi− (k + 1) (51) ⇔ di (k + 1) < (bi (k + 1) − 1) xi− (k + 1) ⇔ ρi (k + 1) > 1 xi (k + 1) < xi− (k + 1) (52) ⇔ di (k + 1) > (bi (k + 1) − 1) xi− (k + 1) ⇔ ρi (k + 1) < 1 xi (k + 1) = xi− (k + 1) (53) ⇔ di (k + 1) = (bi (k + 1) − 1) xi− (k + 1) ⇔ ρi (k + 1) = 1 Appl. Sci. 2021, 11, 9020 20 of 38 for any i ∈ n. We first prove (51) in detail by using contradiction arguments and logics of propositions rules as follows. Let us assume that ρi (k + 1) ≤ 1 for some given i ∈ n. Then, one obtains, from (49), that bi (k + 1) xi− (k + 1) > di (k + 1) ≥ (bi (k + 1) − 1) xi− (k + 1) so that ρi (k + 1) ≤ 1 ⇒ di (k + 1) ≥ (bi (k + 1) − 1) xi− (k + 1) which is the contrapositive logic proposition and then equivalent, logic proposition to the logic implication di (k + 1) < (bi (k + 1) − 1) xi− (k + 1) ⇒ ρi (k + 1) > 1 and the part “⇒” is proved in (51) for any i ∈ n if xi (k + 1) > xi− (k + 1). Now, let us assume that ρi (k + 1) > 1 and di (k + 1) ≥ (bi (k + 1) − 1) xi− (k + 1). Then, ρi (k + 1) = bi (k + 1) − di (k + 1)yi− (k + 1) ≤ bi (k + 1) − (bi (k + 1) − 1) xi− (k + 1)yi− (k + 1) = bi (k + 1) − (bi (k + 1) − 1) xi− (k + 1)yi− (k + 1) = 1 which contradicts that ρi (k + 1) > 1. Then, ρi (k + 1) > 1 ⇒ di (k + 1) < (bi (k + 1) − 1) xi− (k + 1) which proves the part “⇐” of (51) for any i ∈ n. As a result, (51) has been proved for any i ∈ n. The proofs of (52) and (53) are very close to the one above and were therefore omitted. Now, note, from (49), that, for any given i ∈ n, 1 1 =  <1 (54) ρ i ( k + 1) µ i ( k ) µi (k ) bi (k + 1) − di (k + 1)yi− (k + 1) d i ( k +1) µ i ( k ) if µi (k ) bi (k + 1) − di (k + 1)yi− (k + 1) > 1, equivalently, if xi− (k + 1) >  µi (k )bi (k +1)−1 and, in a similar way, one has 1 d i ( k + 1) µ i ( k ) > 1 if xi− (k + 1) < (55) ρ i ( k + 1) µ i ( k ) µ i ( k ) bi ( k + 1 ) − 1 1 d i ( k + 1) µ i ( k ) = 1 if xi− (k + 1) = (56) ρ i ( k + 1) µ i ( k ) µ i ( k ) bi ( k + 1 ) − 1 Extinction at the limit as time tends to infinity needs some extra constraints, as d ( k +1) µ ( k ) follows, since xi (k + 1) < µ (ik)b (k+i1)−1 does not hold at the limit with the left-hand i i side and right-hand side both converging to zero. So, the convergence to extinction of xi (k) as k → ∞ implies that xi (k) = bi (k) xi− (k) − di (k ) → 0 as k → ∞ , then, equivalently, xi−1 (k) = yi (k) → ∞ as k → ∞ and we convey that also xi− (k ) → 0 (basically, the asymp- totic convergences to the right and left of the sampling instants jointly  hold),  bi ( k ) → 1 di ( k ) and di (k) → 0 as k → ∞ with xi ( k ) → 0 as k → ∞ (so that di (k) = o xi−1 (k) as k → ∞ ) d (k) and then ρi (k ) = bi (k ) − xi (k) → 1 as k → ∞ from (49). i The following result, whose proof involves the use of the above considerations, is concerned with non-asymptotic and asymptotic extinction conditions. Theorem 4. The following properties hold: (i) Let us assume that the i − th species has the best fitness, xi (0) > 0 is finite, that the i − th independent consumption fulfills di (k + 1) < bi (k + 1) xi (k + 1); ∀k ∈ N0 and that there is Appl. Sci. 2021, 11, 9020 21 of 38 a strict sequence of non-negative integer numbers {k ` }∞ `=0 ⊂ N0 with k `+1 − k ` ≤ k < + ∞; j=k ∀` ∈ N0 such that ∏ j=k`+1 [ρi (k `+1 + 1 − j)µi (k `+1 − j)] ≥ ρ > 1. Then, there is no ` asymptotic extinction of the i − th species with the best fitness and its solution sequence is bounded for all time. d (k) (ii) Let us assume that (1 − µi )∑nj=1 c ji xi (0) < 1, bi (k ) → 1 , di (k) → 0 and xi (k) → 0 as i k → ∞ and that lim inf µi (k) ≤ 1 with ∑∞ k =1 ( µi ( k ) − 1) = + ∞. Then, there is asymptotic k→∞ extinction of the i-th species with the best fitness and then that of all the remaining ones. Proof. From recursive calculations by using (50), the inverse of the i − th population is given by h i yi (k + 1) = ∏kj=0 (ρi (k + 1 − j)µi (k − j))−1 yi (0) (57) −1 µi ( j)−1 + ∑kj=0 ∏k`= − 1 n K ( j)µ ( j)ρ ( j+1) ∑`=1 `i   j ( i (` ρ + 2 ) i (` µ + 1 )) ( c j j ( ) `i ( )) λ i i i k −1 with the notation convention ∏`= j (ρi (` + 2)µi (` + 1))−1 = 1 for j = k, λii ( j) = 1,   λ`i ( j) ≤ 1; ∀`(6= i ) ∈ n since the i − th species has the best fitness. Since there exists a strictly increasing sequence {k ` }∞`=0 ⊂ N0 with k `+1 − k ` ≤ k < + ∞; ∀` ∈ N0 such that j=k ∏ j=k``+1 [ρi (k `+1 + 1 − j)µi (k `+1 − j)] ≥ ρ > 1 then (58) lim ∏kj=0 ρi (k + 1 − j)µi (k − j)−1   k→∞ and −1 µi ( j)−1 ∑kj=0 ∏k`= −1   j ( ρi (` + 2) µi (` + 1)) Ki ( j ) µ i ( j ) ρ i ( j +1) (59) ×(∑n`=1 c`i ( j)λ`i ( j)) < +∞; ∀k ∈ N0 since xi (0) > 0 is finite and then yi (0) > 0 is finite as well. Then, {yi (k )}∞ k =0 < + ∞ from (57) and { xi (k )}∞ k =0 does not converge to zero. In addition, from that inverse sequence solution formula {yi (k)}∞ k =0 ⊂ R + and then x k ∞ { i ( )}k=0 is bounded since it is convergent. Property (i) has been proved. To prove Property (ii), first note that the joint constraints di (k + 1) < bi (k + 1) xi (k + 1); ∀k ∈ N0 and (1 − µi )∑nj=1 c ji xi (0) < 1, the second one implying also that (1 − µi )∑nj=1 c ji x j (0) < 1 since the i − th species has the best fitness and being always guaranteed if µi ≥; 1, guarantee that { xi (k)}∞ k =0 ⊂ R0+ directly from the time-varying generalization of (2)–(3), even if µi < 1. The constraints bi (k ) → 1 , di (k) → 0 d (k) and xi (k) → 0 as k → ∞ guarantee that eventual asymptotic extinction of the i − th species i jointly holds at the left and right of the sampling instants with ρi (k) → 1 , as k → ∞ , while the constraints, either lim inf µi (k) < 1 or lim inf µi (k) = 1 with ∑∞ k =1 ( µi ( k ) − 1) = + ∞, k→∞ k→∞ ∞ guarantee that {yi (k)}∞ ∞ k =0 → ∞ and then { xi ( k )}k =0 → 0 so that x j ( k ) k =0 → 0 , since  the i − th population has the best fitness; ∀ j ∈ n from the formula which obtains the solution of the inverse population from the initial conditions.  3. Stability Results Obtained by Lyapunov Method It is now assumed that the Beverton–Holt equation is time-varying, in general, rather than parameterized by constant parameters. A Lyapunov sequence candidate is proposed for stability analysis around eventual non-extinction equilibrium points. From Assumption d (k) 2, an artificial auxiliary sequence {νi (k)}k∈N ⊂ R is defined by νi (k) = x−i (k) ∈ (−∞, bi (k)]; i Appl. Sci. 2021, 11, 9020 22 of 38 ∀i ∈ n, ∀k ∈ N (note that di (k) is not defined for k = 0). That definition, with the above constraint, guarantees that di (k + 1) = νi (k + 1) xi− (k + 1)≤ bi (k + 1) xi− (k + 1); ∀i ∈ n, ∀k ∈ N (60) Remark 3. The above artificial sequence allows an easy integration of the eventual independent consumption contribution to the stability of non-extinction equilibrium points. The above artificial sequence is not adequate to describe finite-time extinction xi (k) = xi− (k + 1) = 0 with di (k + 1) < 0 (individuals of the species i − th enter the habitat from outside) implying that the extinguished i − th species to the left of the (k + 1) sampling instants recovers from previous extinction to the right of the (k + 1) sampling instant.  − We ∞ can now describe the model to the right of the sampling instants as follows, if xi (k) k=0 ⊂ R+ , ∀i ∈ n, µi (k)Ki (k )(bi (k + 1) − νi (k + 1)) x i ( k + 1) = x (k); ∀i ∈ n, ∀k ∈ N (61) Ki (k ) + (µi (k) − 1)∑nj=1 c ji (k ) x j (k ) i di ( k ) for given initial finite conditions xi (0) > 0; ∀i ∈ n, where νi (k) = xi− (k ) ∈ (−∞, bi (k)]. Let us define the i − th population inverse yi (k) = xi−1 (k); ∀i ∈ n, ∀k ∈ N so that yi (k + 1) = (µi (k )(bi (k + 1) − νi (k + 1)))−1 yi (k)   (62) + µ (k)K (k)(µbi ((kk)− 1 n +1)−ν (k+1)) ∑ j=1 ji ( ) j ( ) i ( ) c k x k y k i i i i which can be rewritten as    µ (k )−1 yi (k + 1) = ( ai (k) + hi (k))yi (k ) = ai (k) 1 + iK (k) ∑nj=1 c ji (k ) x j (k ) yi (k) i    (63) µi (k )−1 µi (k )−1 = ai (k) 1 + K (k) + K (k) ∑nj(6=i)=1 c ji (k) x j (k) yi (k) i i where ai (k) = (µi (k )(bi (k + 1) − νi (k + 1)))−1 (64)   µ (k )−1 hi (k) = µ (k)K (k)(bi (k+1)−ν (k+1)) ∑nj=1 c ji (k) x j (k) i i i i   (65) (µi (k)−1) ai (k) n = K (k) i ∑ j = 1 c ji ( k ) x j ( k ) The following result relies on the asymptotic stability of the non-extinction equilibrium point of the competition impulsive Beverton–Holt equation. Theorem 5. Let us assume that all the competing species have non-zero equilibrium points xi ∈ R+ . Those equilibrium points are asymptotically stable, provided that the following constraints hold:    (a) ai (k) + hi (k) = ai (k) 1 + µiK(k()− k) 1 ∑ n j = 1 c ji ( k ) x j ( k ) < 1; i ∀ k ∈ N0 h   i (b) lim sup ai (k) 1 + µiK(k()− k ) 1 ∑ n j = 1 c ji (k) x j (k) ≤1 i k→∞ Appl. Sci. 2021, 11, 9020 23 of 38 and that, furthermore, at least one of the subsequent four constraints hold for some set of real parameters θi > 0; ∀i ∈ n: (c) xi (k) < (1 + ai (k) + hi (k)) x2i ; ∀k ∈ N0 (d) ∑in=1 θi xi ( k ) (1 + ai (k) + hi (k)) > 2∑in=1 θi xi ; ∀ k ∈ N0 max xi2 (k ) 2 1≤ i ≤ n ∑in=1 θi (1−( ai (k )+hi (k )) ) (e) min xi (k ) < 2∑in=1 θi [1− ai (k )−hi (k )] xi 1≤ i ≤ n min (1−( ai (k )+hi (k ))2 ) ∑in=1 θi xi2 (k ) 1≤ i ≤ n (f) ∑in=1 θi xi (k ) < 2 max (1− ai (k )−hi (k )) xi 1≤ i ≤ n Proof. Let us define yi = xi−1 ; ∀i ∈ n with xi being a nonzero equilibrium point of the i-th species and yei (k) = yi (k) − yi ; ∀i ∈ n, ∀k ∈ N0 . Let us consider the Lyapunov sequence candidate V (k) = ∑in=1 θi (yi (k) − yi )2 ; ∀k ∈ N0 for some set of design scalar coefficients θi > 0; ∀i ∈ n. Then, the one-step incremental Lyapunov sequence candidate verifies the subsequent chain of relations: ∆V (k) = V (k + 1) − V (k) = ∑in=1 θi ye2i (k + 1) − ye2i (k)  = ∑in=1 θi (yei (k + 1) + yei (k)) (yei (k + 1) − yei (k)) = ∑in=1 θi ([ ai (k) + hi (k)]yi (k) − yi + yi (k) − yi ) ×([ ai (k) + hi (k)]yi (k) − yi − yi (k) + yi ) = ∑in=1 θi ([ ai (k) + hi (k) + 1]yi (k) − 2yi ) ([ ai (k) + hi (k) − 1]yi (k)) (66) = ∑in=1 θi [ ai (k) + hi (k) + 1] [ ai (k) + hi (k) − 1]y2i (k) −2∑in=1 θi yi [ ai (k) + hi (k) − 1]yi (k) = −∑in=1 θi ((1 + ai (k) + hi (k)) yi (k) − 2yi )(1 − ai (k) − hi (k))yi (k) ≤ −∑in=1 ((1 + ai (k) + hi (k)) yi (k) − 2yi ) min [θi (1 − ai (k) − hi (k))yi (k)]; 1≤ i ≤ n ∀ k ∈ N0 Other upper-bounds of ∆V (k) are clearly the following ones:  h i − ∑in=1 θi 1 − ( ai (k) + hi (k))2 min y2i (k) 1≤ i ≤ n (67) +2(∑in=1 θi yi [1 − ai (k) − hi (k)]) max yi (k) 1≤ i ≤ n and − min 1 − ( ai (k) + hi (k))2 ∑in=1 θi y2i (k)   1≤ i ≤ n (68) +2 max (1 − ai (k) − hi (k))(∑in=1 θi yi (k)yi ) 1≤ i ≤ n Thus, if the constraints (a) and (b) and one of the constraints from (c) to (f) hold, then {−∆Vk }∞ ∞ ∞ k =0 ⊂ R+ , so that {Vk }k =0 is bounded and strictly decreasing, so that {Vk }k =0 → 0 −1 so that yi (k) → yi = xi , equivalently, xi (k) → xi ; ∀i ∈ n; ∀k ∈ N0 .  Remark 4. Note that the constraint (a) of Theorem 5 requires, as a necessary condition, µi (k)(bi (k + 1) − νi (k + 1)) > 1; ∀i ∈ n; ∀k ∈ N0+ . Appl. Sci. 2021, 11, 9020 24 of 38 Note also that a sufficient condition for Theorem 5 to hold is ∑in=1 θi ((1 + ai (k) + hi (k )) yi (k) − 2yi )(1 − ai (k) − hi (k )) > 0; (69) ∀i ∈ n, ∀k ∈ N0 although it is not as easy to test as the conditions given in the theorem statement. Note that stability around the zero equilibrium implies the asymptotic extinction of all the species, This also implies that the inverse populations diverge asymptotically to +∞. Note that 1 y i ( k + 1) = µi (k )Ki (k )(bi (k +1)−di (k +1)yi− (k +1)) ( Ki ( k ) y i ( k ) + µ i ( k ) − 1) µi (k )−1  y (k)  (70) + µi (k )Ki (k )(bi (k +1)−di (k+1)yi− (k +1)) ∑nj(6=i)=1 c ji (k) yij (k) The asymptotic extinction is formalized through the subsequent result. Theorem 6. Let us assume that the following conditions hold: (a) {Ki (k)}∞ ∞ h i ∞  k =0 ⊂ 0 , K , { µi ( k )}k =0 ⊂ (1, µ ), { bi ( k )}k =0 ⊂ 0, b , cii ( k ) = 1, ∞ c ji (k ) k=0 ⊂ [0, 1] and di (k) ∈ −∞, bi (k) xi− (k) ; ∀i ∈ n, ∀k ∈ N0 .   (b) If there  k ∈ N0 such that xi (k) = 0, then xi (k + 1) = −di (k + 1) = 0  is some finite (c) lim ∏kj=0 [αi ( j)] = +∞ k→∞ Then, { xi (k)}∞ k=0 → 0 ; ∀i ∈ n if Ki ( k ) y i ( k ) + µ i ( k ) − 1 lim inf  >1 k→∞ µi (k)Ki (k) bi (k + 1) − di (k + 1)yi− (k + 1) Proof. Note that the assumption di (k ) ≤ bi (k ) xi− (k) is equivalent to bi (k) − di (k )yi− (k) ≥ 0; ∀i ∈ n, ∀k ∈ N0+ and implies that xi (k) ≥ 0 and yi (k) ≥ 0; ∀i ∈ n, ∀k ∈ N0 . Let us define the auxiliary variables for i ∈ n, k ∈ N0 as Ki ( k ) y i ( k ) + µ i ( k ) − 1 αi ( k ) = (71) µi (k )Ki (k) bi (k + 1) − di (k + 1)yi− (k + 1)  µi ( k ) − 1 β i (k) = (72) µi (k )Ki (k) bi (k + 1) − di (k + 1)yi− (k + 1)  n y (k) γi ( k ) = ∑ j(6=i)=1 c ji (k) yij (k) (73) so that yi (k + 1) = αi (k )yi (k) + β i (k)γi (k); ∀i ∈ n, ∀k ∈ N0 (74) One obtains, via recursive calculations with the above relation, that     yi (k + 1) − ∑kj=0 ∏kj=0 [αi (`)]k− j β i ( j)γi ( j) = ∏kj=0 [αi ( j)] yi (0); (75) ∀i ∈ n, ∀k ∈ N0 And, if yi (0) > 0 and finite (equivalently, 0 < xi (0) < +∞) for all i ∈ n, then    k k lim yi (k + 1) − ∑ j=0 ∏ j=0 [αi (`)]k− j β i ( j)γi ( j) = +∞; ∀i ∈ n (76) k→∞ The following cases can arise: Appl. Sci. 2021, 11, 9020 25 of 38 Case 1: {yi (k )}∞ ∞ k =0 → + ∞ ; ∀i ∈ n, equivalently, { xi ( k )}k =0 → 0 ; ∀i ∈ n. The proof follows directly. Case 2: {yi (k)}∞    k k k− j β ( j)γ ( j) = +∞ for k =0 is bounded and lim ∑ j =0 ∏ j =0 [ α i (`)] i i k→∞ some i ∈ n. Case 2 can hold through several sub-cases as follows: a) Sub-case 2a: Either {αi (k)}∞ ∞ k=0 → + ∞ or { β i ( k )}k =0 → + ∞ . Since the numerators of their respective general terms are finite, then either the denominators converge to zero or they have some element which is zero, that is, either bi (k) − di (k)yi− (k )(≥ 0) → 0 as k → ∞ for some i ∈ n or there are some k ∈ N0 and some i ∈ n, such that bi (k) − di (k)yi− (k) = 0. However, these two sub-cases imply, respectively, { xi (k )}∞ k =0 → 0 and xi (k) = 0, which implies that xi− (k + 1) = 0 and xi (k + 1) = −di (k + 1) = 0 from the assumption (b). In both cases, extinction is achieved either asymptotically or in a finite time and the claim of Case 2 that {yi (k)}∞ k=0 is bounded cannot hold for any i ∈ n so that {yi (k)}∞ k =0 → + ∞ ; ∀ i ∈ n. Thus, Case 2, supported by Sub-case 2a, cannot hold. ∞ b) Sub-case 2b: Either {γi (k)}∞ k =0 → + ∞ or γ j ( k ) k =0 → + ∞ for some given i ∈ n ∞ y (k) ∞ n o → ∞,  and some j(6= i ) ∈ n. Since the sequences c ji (k) k=0 ⊂ [0, 1] then either yi (k) j k =0 o∞ = +∞, for j(6= i ) ∈ n. Since {yi (k)}∞ n y (k) or yi (k) k =0 is claimed to be bounded, then the j k =0 denominators of the quotients defining the general terms of the above sequence either converge to zero as k → ∞ or are zero for some finite k ∈ N0 . However, this would imply the unboundedness of some of the subpopulations contradicting the boundedness property of Theorem 2, which still holds for the extended Beverton–Holt evolving under the time-varying parameterizing sequences {Ki (k )}∞ ∞ ∞ ∞ k =0 , { µi ( k )}k =0 , { bi ( k )}k =0 and { di ( k )}k=0 , fulfilling the assumption (a) of this theorem. Therefore, Case 2 cannot hold when being supported by Sub-case 2b.    c) Sub-case 2c: lim ∑kj=0 ∏kj=0 [αi (`)]k− j β i ( j)γi ( j) = +∞ with k→∞ = +∞ and { β i ( j)γi ( j)}∞    lim ∑kj=0 ∏kj=0 [αi (`)]k− j j=0 being bounded for any i ∈ n. k→∞   Since lim ∏kj=0 [αi ( j)] = +∞ and since the above partial sum diverges to +∞ as k → ∞ , k→∞   one also concludes, from the assumed condition (c), i.e., lim ∏kj=0 [αi ( j)] = +∞, that k→∞ the Sub-case 2c reduces to Case 1. As a conclusion, only Case 1, including Sub-Case 2c, is possible, which implies the asymptotic extinction of all the subspecies.  Remark 5. As it could be expected from intuition, a large hunting/fishing quota or a certain level of positive independent consumption might lead to extinction. Such conclusions arise from inspecting either the zeroing in finite time or the asymptotic convergence to zero of the amount bi (k + 1) − di (k + 1)yi− (k + 1). In a parallel way, sufficiently small hunting/fishing quotas or small positive independent consumption can keep the populations free of extinction. 4. About Extinction Conditions under Small Intrinsic Growth Rates Note that Theorem 6 gives extinction conditions for intrinsic growth rates exceeding unity, which is an usual condition in the background literature on Beverton–Holt equa- tions [5–12,18,26–29,33]. In this section, the asymptotic extinction is studied for the case of intrinsic growth rates being less than unity. Remark 6. Note that extinction at the right of the k − th sampling instant of the i − th species oc- curs if di (k ) = bi (k ) xi− (k) and it occurs asymptotically if di (k) − bi (k ) xi− (k) → 0 as k → ∞ . This can happen irrespective if the intrinsic growth rate exceeds unity or is less than unity. The first case, which is the typical situation in the Beverton–Holt equation, might be addressed as a particular case included in Theorem 6. The condition is difficult to achieve in practical cases, since it requires an exact monitoring of the independent consumption through time. In the following and for the case of intrinsic growth rates below unity, this case is excluded from the analysis. Appl. Sci. 2021, 11, 9020 26 of 38 To fix ideas, we first discuss the case when just a single species is involved. For this purpose, we establish two preliminary necessary technical results for possible extinction either in finite time or asymptotic extinction, while keeping the solution sequence non- negative and bounded for all time under non-negative initial conditions. Proposition 5. Let us consider one single-species impulsive Beverton–Holt equation, Equation (1). Let us assume that {µ(k)}∞ ∞ ∞ k =0 ⊂ (0, 1], { b ( k )}k =1 ⊂ (0, 1] and { d ( k )}k =1 ⊂ (− ∞, + ∞ ) Then, for any given k ∈ N0 , x (k) = 0 ⇒ x − (k + 1) = 0 and, furthermore, x (k + 1) = 0 if and only if d(k + 1) = 0. If { x (k)}∞ ∞ ∞ k =0 → 0 , then { d ( k )}k =1 → 0 and, if, in addition, { b ( k )}k =1 ⊂ R+ , − ∞ then { x (k)}k=0 → 0 . Proof. One obtains, from (1), that [b(k+1)µ(k)K (k)+d(k+1) (1−µ(k))] x (k)−d(k+1)K (k) x ( k + 1) = K (k )−(1−µ(k )) x (k ) ; (77) ∀ k ∈ N0 Since K (k) > 0, x (k) = 0 ⇒ x − (k + 1) = 0 from (1) and x (k + 1) = 0 if and only x (k +1)+d(k +1) if d(k + 1) = 0 and, if b(k + 1) > 0, then x − (k + 1) = b ( k +1) = 0. Note that, if { x (k)}∞ k =0 → 0 , then b ( k +1) µ ( k ) K ( k ) x ( k ) d(k + 1) = K (k)−(1−µ(k)) x(k) − x (k + 1) → 0 as k→∞ and also {b(k)}∞ { x (k)}∞ x (k +1)+d(k +1) x − ( k + 1) = b ( k +1) → 0 as k → ∞ if k =1 ⊂ R+ since k =0 → 0 and {d(k)}∞ k =1 → 0.  Proposition 6. Let us consider a Beverton–Holt equation for one single species. Let us assume that {µ(k )}∞ ∞ ∞ k =0 ⊂ (0, 1], { b ( k )}k =1 ⊂ (0, 1] and { d ( k )}k =1 ⊂ (− ∞, + ∞ ). Then, for any given k ∈ N0 , x (k + 1) ∈ [0, +∞), with x (k) ≥ 0, if and only if h    d ( k +1) K ( k ) K (k) x (k) ∈ min 0, b(k+1)µ(k)K (k)+d(k+1)(1−µ(k)) , 1−µ(k) ; (78) ∀ k ∈ N0 Proof. Let us consider several cases in h (77) with the given assumptions, namely,  d ( k +1) K ( k ) K (k) Case a: x (k + 1) ≥ 0, if x (k) ∈ , b(k +1)µ(k )K (k )+d(k +1)(1−µ(k )) 1−µ(k ) . Note that the limits of the admissible domain of x (k) are mutually compatible, since d ( k +1) 1 b(k +1)µ(k )K (k )+d(k +1)(1−µ(k )) < 1− µ ( k ) b ( k +1) µ ( k ) K ( k ) ⇔ d ( k + 1) < d ( k + 1) + 1− µ ( k ) b ( k +1) µ ( k ) K ( k ) which trivially holds, since 1− µ ( k ) > 0. Since x (k) is also requested to be non- negative, then (78) is used. It was proved that (78) implies that x (k + 1) ≥ 0. Now, note that the conditions {µ(k)}∞ ∞ ∞ k =0 ⊂ (0, 1], { b ( k )}k =1 ⊂ (0, 1] and { d ( k )}k =1 ⊂ (− ∞, + ∞ ) also imply that the domain of x (k ) is bounded and, from (77), x (k + 1) is also bounded. Thus, Case a fully verifies the statement  claim.  K (k) d ( k +1) K ( k ) Case b: x (k + 1) ≥ 0, if x (k) ∈ , 1−µ(k ) b(k +1)µ(k )K (k )+d(k +1)(1−µ(k )) . However, a necessary condition for the above domain of values of x (k) to be well-posed is that 1 d ( k +1) 1− µ ( k ) < b(k +1)µ(k )K (k )+d(k +1)(1−µ(k )) b ( k +1) µ ( k ) K ( k ) ⇔ d ( k + 1) > d ( k + 1) + 1− µ ( k ) Appl. Sci. 2021, 11, 9020 27 of 38 which is a contradiction. Then, Case b does not verify the statement claim. K (k) d ( k +1) K ( k ) Case c: x (k + 1) ≥ 0, if x (k ) = 1−µ(k) = b(k+1)µ(k)K (k)+d(k+1)(1−µ(k)) and the indeter- mination 0/0, resulting from (77), is solvable with a positive limit, as the numerator and denominator converge to zero. However, the identity of both confluent values of x (k) leads to the following contradiction: b ( k + 1) µ ( k ) K ( k ) 0 = d ( k + 1) − d ( k + 1) = >0 1 − µ(k) so that Case c does not verify the statement claim. As a result, only Case a satisfies the statement claim and the proof is complete.  The next main result of this section concerned with the case in which one single species relies on the asymptotic convergence of the solution to extinction when the solution to the right of the sampling instants is a non-negative strictly decreasing sequence, that is, { x (k) − x (k + 1)}∞ ∞ k =0 ⊂ R+ and { x ( k )}k =0 (⊂ R+ ) → 0 . Theorem 7. Let us consider the Beverton–Holt equation for one single species. Let us assume that, if x (k0 ) = 0 for any finite k0 ∈ N0 , then d(k0 + 1) = 0. Let us also assume that 0 < x (0) = x − (0) < 1−µ(0) , {µ(k)}∞ ∞ ∞ K (0) k =0 ⊂ (0, 1], { K ( k )}k =1 ⊂ 0, K with K < + ∞, { b ( k )}k =1 ⊂  (0, 1]. Let us also assume that  i [(b(k+1)µ(k)−1)K (k)+(1−µ(k)) x (k)] x (k) b ( k +1) µ ( k ) K ( k ) x ( k ) d ( k + 1) ∈ K (k )−(1−µ(k )) x (k ) , K (k )−(1−µ(k )) x (k ) ; (79) ∀ k ∈ N0 ∞ Then, { x (k)}∞ − ∞ k =0 ⊂ R0+ , { x ( k )}k =0 ⊂ R0+ , { d ( k )}k =1 ⊂ R are bounded, with ∞ { x − (k)}k=0 → 0 and { x (k)}k=0 → 0 being strictly decreasing and {d(k)}∞ ∞ k =1 → 0 . Proof. Let us consider the Lyapunov sequence candidate V (k) = x (k ); ∀k ∈ N0 so that K (k) ∆V (k) = V (k + 1) − V (k) = x (k + 1) − x (k ). Then, from Proposition 6, 0 < x (k) < 1−µ(k) for any k ∈ N0 implies that x (k + 1) ≥ 0, so that (K (k) − (1 − µ(k)) x (k))∆V (k) = (K (k) − (1 − µ(k)) x (k)) (V (k + 1) − V (k)) = [(b(k + 1)µ(k) − 1)K (k) + d(k + 1) (1 − µ(k)) + (1 − µ(k)) x (k)] x (k) − d(k + 1)K (k) (80) = [(b(k + 1)µ(k) − 1)K (k) + (1 − µ(k)) x (k)] x (k) − d(k + 1) [K (k) − (1 − µ(k)) x (k)]≤ 0; ∀ k ∈ N0 and { x (k)}∞ K (0) k =0 ⊂ R0+ is bounded and non-increasing, since 0 < x (0) < 1− µ (0) if [(b(k + 1)µ(k) − 1)K (k) + (1 − µ(k)) x (k)] x (k) d ( k + 1) ≥ ; ∀ k ∈ N0 (81) K (k ) − (1 − µ(k)) x (k ) which implies that ∆V (k) ≤ 0; ∀k ∈ N0 . Since K (k ) > (1 − µ(k)) x (k ), it follows, from (80), that x (k + 1) < x (k), so that { x (k)}∞ k =0 is strictly decreasing, if x ( k ) > 0 and (81) is a strict inequality, ∆V (k) → 0 as k → ∞ and that x (k) → 0 as k → ∞ . If x (k + 1) ≥ 0, then b ( k + 1) µ ( k ) K ( k ) x ( k ) d ( k + 1) ≤ b ( k + 1) x − ( k + 1) = ; ∀k ∈ N0 (82) K (k) − (1 − µ(k)) x (k) ∞ and { x − (k)}k=0 ⊂ R0+ is also bounded, since x (0) = x − (0) is finite and {d(k )}∞ k =1 ⊂ (−∞, +∞). Combining both inequalities, (81), in its strict version, and (82), yields (79). This Appl. Sci. 2021, 11, 9020 28 of 38 also implies that x (k) → 0 as k → ∞ so that d(k + 1) → 0 and x − (k + 1) → 0 as k → ∞ .  Remark 7. Note that Theorem 7 holds, in particular and as stated, if b(k + 1) = µ(k) ≡ 0 since it involves, for all k ∈ N0 the condition − x (k) < d(k + 1) ≤ 0, implying that x (k + 1) − x (k) = −d(k + 1) − x (k) < x (k) − x (k) = 0 if x (k) 6= 0, so that the one-step increment of the Lyapunov sequence is negative. Let us consider now the extended competition model from (4) for the extinction of some species i ∈ n in the competition model with possible time-varying parameterization to conclude that x i ( k +1) xi ( k )    d i ( k +1) (83) = [bi (k + 1)µi (k)Ki (k) + di (k + 1)(1 − µi (k))] − xi ( k ) Ki (k ) + ∑nj(6=i)=1 c ji (k )(µi (k ) − 1) x j (k ) ∀ k ∈ N0 1 × < 1; ∀k ∈ N0 Ki (k ) − (1 − µi (k )) xi (k ) + ∑nj(6=i)=1 c ji (k)(µi (k) − 1) x j (k) with the positive denominator in (83) implying a non-negative numerator, implying, equivalently, that 0 ≤ xi (k + 1) < xi (k) if xi (k) > 0 for all k ∈ N0 if the “ad-hoc” version of (79) holds in the form bi ( k + 1 ) µ i ( k ) K i ( k ) x i ( k ) d i ( k + 1) ≤ Ki (k )−(1−µi (k )) xi (k )+∑nj(6=i)=1 c ji (k )(µi (k )−1) x j (k ) ; (84) ∀ k ∈ N0 (bi (k+1)µi (k)−1)Ki (k)+(1−µi (k)) xi (k)−∑nj(6=i)=1 c ji (k)(µi (k)−1) x j (k) d i ( k + 1) > Ki (k )−(1−µi (k )) xi (k )+∑nj(6=i)=1 c ji (k )(µi (k )−1) x j (k ) x i ( k ); (85) ∀k ∈ N0 which leads to the subsequent result. Theorem 8. Let us consider a competition Beverton–Holt equation for n species. Let us assume that, if xi (k0 ) = 0 for any finite k0 ∈ N0 and i ∈ n then di (k0 + 1) = 0. Let us also assume that x j (0) = x − j (0) ≥ 0; ∀ j (6 = i ) ∈ n; Ki (0)+∑nj(6=i)=1 c ji (0)(µi (0)−1) x j (0) 0 < xi (0) = xi− (0) < 1− µ i (0) >, ∞ {µi (k)}∞ ∞ (1, µ` ) ⊂ (1, +∞), b j (k)  k =0 ⊂ (0, 1], { µ` ( k )}k =0 ⊂ k =1 ⊂ (0, 1], ∞ K j (k) k=1 ⊂ 0, K j ⊂ (0, +∞) and −∞, b j (k) x −   d j (k) ⊂ j (k) ]; ∀ j, `(6= i ) ∈ n; ∀k ∈ N. Let us also assume that (bi (k+1)µi (k)−1)Ki (k)+(1−µi (k)) xi (k)−∑nj(6=i)=1 c ji (k)(µi (k)−1) x j (k)  d i ( k + 1) ∈ Ki (k )−(1−µi (k )) xi (k )+∑nj(6=i)=1 c ji (k )(µi (k )−1) x j (k ) x i ( k ),  (86) bi ( k + 1 ) µ i ( k ) K i ( k ) x i ( k ) K (k )−(1−µ (k )) x (k )+∑n c (k )(µ (k )−1) x (k ) ; ∀k ∈ N0 i i i j(6=i )=1 ji i j ∞ Then, { xi (k)}∞ ⊂ R0+ , xi− (k) k=0 ⊂ R0+ , {di (k)}∞  k = 0 k =1 ⊂ R are bounded, with ∞ xi (k ) k=0 → 0 and { xi (k )}k=0 → 0 being strictly decreasing and {di (k)}∞ ∞  − k =1 → 0 . {x (k )} − i ∞ k =0 → 0 and {xi (k )}∞k = 0 → 0 being strictly decreasing and {d i (k )}k =1 → 0 . ∞ Appl. Sci. 2021, 11, 9020 5. Numerical Experiments 29 of 38 This Section aims at illustrating, through numerical simulation examples, some of the theoretical results discussed in the previous Sections. In this way, this Section is organized into a 5. number of examples Numerical concerned with the Theorems and results stated in Sections 2– Experiments 4. This Section aims at illustrating, through numerical simulation examples, some of the theoretical results discussed in the previous Sections. In this way, this Section is organized Example 2. Let us consider the system (1) parameterized by b(k) = 1, d(k) = 0, K(k) = 300 into a number of examples concerned with the Theorems and results stated in Sections 2–4. and μ( k ) = 0.6 + 5 × 0.6 k with x(0) = 100. The values of μ (k ) with iteration along with the evolution Exampleof2.x(k) areconsider Let us displayed in Figure the system (1) 1. Thus, Lemma parameterized 1 holds, by b(k) since = 1, d(k) = the growth 0, K(k) = 300 and rate converges µ(k ) = 0.6to+0.6 0.6k μ 5 ×and with ∈ ( 0= ( k )x(0) ,1 )100. forThe k ≥ 5 .ofItµis all values (kobserved, in Figure ) with iteration 1, how along with the the evolution of x(k) extinguishes population are displayed ininFigure finite 1. Thus, time, Lemma since x(j) 1= holds, 0 for j>since 30 the growth rate according converges to Lemma 1. to On 0.6 and ) for allμk( ≥ k the other ∈ (0, 1when µ(k ) hand, = 1+ k ) 5. It is observed, 0.62 , we have that1, μ in Figure how (k) the population converges to extinguishes unity (in fact in finite ( ) = + k μ (k) is forced to be exactly unity from j = 10 onwards) and we are in the condition of , we time, since x(j) = 0 for j> 30 according to Lemma 1. On the other hand, when µ k 1 0.62 have that µ(k) converges to unity (in fact µ(k) is forced to be exactly unity from j = 10 onwards) applying Lemma 2. Therefore, Figure 2 shows how the population remains strictly posi- and we are in the condition of applying Lemma 2. Therefore, Figure 2 shows how the population tive for all discrete time (and even converges to a positive constant value) as Lemma 2 remains strictly positive for all discrete time (and even converges to a positive constant value) as claims. Lemma 2 claims. 6 (k) 5 4 3 2 1 0 0 5 10 15 20 25 30 35 Iteration 300 x(k) 250 200 150 100 50 0 0 5 10 15 20 25 30 35 Iteration FigureFigure 1. Evolution of μ (ofk µ(k) 1. Evolution ) and x(k)with andx(k) with iteration iteration in in Example Example22when μ (converges whenµ(k) k ) converges to a limit to abe-low unity. limit below unity. Example 3. Let us consider the system (1) parameterized by the constant values b(k) = 0.9, d(k) = 0.1, K(k) = 300 and µ(k) = 2.4 with x(0) = 100. The evolution of x(k) is displayed in Figure 3. It can be seen, in this figure, that the population converges to the equilibrium point given numerically by 248.385. This value is in accordance with the calculations in Theorem 1(i), since µ(k ) = 2.4 > max(1, 1/b(k)) = 1.11 and d = 0.1 satisfies the constraint " √ 2 # " √ 2 ! µb−1 K µb+1 K 0.1 = d ∈ − x1 , µ −1 ∪ µ −1 , ∞ = [−0.0863, 47.2741] ∪ [1307, ∞). Appl. Sci. 2021, 11, x FOR PEER REVIEW 35 of 44 Appl. Sci. 2021, 11, 9020 30 of 38 2 (k) 1.9 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0 5 10 15 20 25 30 35 Iteration 170 160 x(k) 150 140 130 Appl. Sci. 2021, 11, x FOR PEER REVIEW 120 36 of 44 110 (μb − 1)K = 300 , as it can be readily seen in Figure 4. In addition, when b = μ −1 (and x2 =100 0 5 10 15 20 25 30 35 μ −1 Iteration d = 0), the 2. trajectory of the population isiteration depicted in Figure 5. It can converges be observed, in this Figure Figure figure, that Evolution 2. Evolution the of μ ( of k )µ(k) population andendsx(k)x(k) and with up with iteration in Example in Example extinguishing, as 2 whenμµ(k) 2 when predicted by ( Theorem toto k ) converges 1(ii). unity unityin finite time. in finite time. Example 280 3. Let us consider the system (1) parameterized by the constant values b(k) = 0.9, 260 K(k) = 300 and μ ( k ) = 2.4 with x(0) = 100. The evolution of x(k) is displayed in d(k) = 0.1, Figure 3. It can be seen, in this figure, that the population converges to the equilibrium 240 point given numerically by 248.385. This value is inx(k) accordance with the calculations in - x (k) Theorem 220 1(i), since μ ( k ) = 2.4 > max( 1, 1 / b( k )) = 1.11 and d = 0.1 satisfies the constraint 200  0.1 = d ∈ − x1 , 2 ∪ ( μb − 1 K   μb + 1 K  2 ), ∞ ( ) 180  μ −1   μ −1      160 = [− 0.0863, 47.2741] ∪ [1307 , ∞ ) . Thus, 140 Theorem 1(i) holds. Furthermore, Theorem 1(i) provides two equilibrium given by x1 = 0.0863 and x2 = 248.385 . Moreover, Figure 3 also depicts the left points 120 value ofthe population x−( k ) . It satisfies at the equilibrium 100 0 5 x + d10 248.385 + 0.1 15 20 25 30 35 276.09 = x − ( k ) = 2 = Iteration = 276.09 , as Theorem 1(i) states. If we now con- b 0.6 sider Figure the Figure 3. Evolution 3.case whenofb(k) Evolution of=x(k) x(k) for for constant 1 constant and parameters d = parameters 0, i.e., there in isinno Example Example 3, including harvesting harvesting norharvesting 3, including and and independent external con- external consumption. consumption. The The population population converges converges to a to a positive positive equilibrium equilibrium point. point. sumption, then the left and right values of the population are the same, as Figure 4 shows, while the equilibrium point is in accordance with Theorem 1(ii), since 300 x(k) 280 - x (k) 260 240 220 280 260 240 x(k) - Appl. Sci. 2021, 11, 9020 x (k) 31 of 38 220 200 180 Thus, Theorem 1(i) holds. Furthermore, Theorem 1(i) provides two equilibrium points given 160 by x1 = 0.0863 and x2 = 248.385. Moreover, Figure 3 also depicts the left value of the population x − (k). It satisfies at the equilibrium 276.09 = x − (k) = x2b+d = 248.385 0.6 +0.1 = 140 276.09, as Theorem 1(i) states. If we now consider the case when b(k) = 1 and d = 0, i.e., there 120 is no harvesting nor independent consumption, then the left and right values of the population 100 are the same, as Figure 4 shows, while the equilibrium point is in accordance 0 5 10 15 20 ( µb −25 1) K 30 35 with Theorem 1(ii), since x2 = µ−1 = 300, as it can be readily seen in Figure 4. In Iteration addition, when b = µ−1 (and d = 0), the trajectory of the population is depicted in Figure 5. Figure 3. Evolution of x(k) for constant parameters in Example 3, including harvesting and external It can be observed, in this figure, that the population ends up extinguishing, as predicted consumption. The population converges to a positive equilibrium point. by Theorem 1(ii). 300 x(k) 280 - x (k) 260 240 220 200 180 160 140 120 100 0 5 10 15 20 25 30 35 Iteration Appl. Sci. 2021, 11, x FOR PEER REVIEW 37 of 44 Figure Figure 4. Evolution 4. Evolution of x(k) of x(k) for constant for constant parameters parameters in Example in Example 3, while 3, while no no harvesting harvesting nornor independent independ- consumptionisisconsidered. ent consumption considered. 180 x(k) 160 x - (k) 140 120 100 80 60 40 20 0 0 20 40 60 80 100 120 140 160 180 200 Iteration Figure 5. Evolution of x(k) for constant parameters in Example 3 and b = −1 with d = 0. −1 µ Figure 5. Evolution of x(k) for constant parameters in Example 3 and b = μ with d = 0. Example Example Letconsider 4. us 4. Let us considerthethe multi-speciesBeverton–Holt multi-species Beverton–Holt Equations Equations(2) (2)and and(3) (3)parameterized param- eterized by the constant values bi (k ) = 1 , d i (k ) = 0 , K1 (k ) = 300 , K 2 (k ) = 200= and by the constant values b i ( k ) = 1, d i ( k ) = 0, K 1 ( k ) = 300, K 2 ( k ) = 200 and K 3 ( k ) 100 and µ1 (k) = 2.4, µ2 (k) = 2 and µ3 (k) = 1.8 with x1 (0) = 100, x2 (0) = 50 and x3 (0) = 30. The K 3 (k ) = 100 and μ1 (k ) = 2.4 , μ 2 (k ) = 2 and μ3 (k ) = 1.8 with x1 (0 ) = 100 , x2 (0) = 50 coupling matrix is given by and x3 (0) = 30 . The coupling matrix is given  by 1 0.4 0.3  C =  0.25 1 0.4   1 0.4 0.3   0.5 0.35  1 C =  0.25 1 0.4   0.5 0.35 1    These parameters correspond to the case when the species 1 dominates the other two. Thus, we are in the condition of applying Propositions 1 and 2 and Theorem 2(i) to guar- antee that all populations are non-negative and, furthermore, x1 (k ) > x2 (k ) > x3 (k ) , as Appl. Sci. 2021, 11, 9020 32 of 38 These parameters correspond to the case when the species 1 dominates the other two. Thus, we are in the condition of applying Propositions 1 and 2 and Theorem 2(i) to guarantee that all populations are non-negative and, furthermore, x1 (k) > x2 (k ) > x3 (k ), as Proposition 2 ensures. These results can be observed in Figure 6, where the evolution of the three populations is displayed. In addition, Figure 7 displays the evolution of the three species when we consider the same parameters as before, but b1 = 0.98, b2 = 0.96 and b3 = 0.95 with d1 = 0.01, d2 = 0.02 and d3 = 0.03 (i.e., b1 > b2 > b3 while d3 > d2 > d1 ). It is checked, in this case, that both Propositions 1 and 2 hold when a harvesting quota, along with an external consumption, is included in the system. It can be also observed, in Figure 7, that species 3 extinguishes. This situation is also predicted by Theorem 2(ii), since Appl. Sci. 2021, 11, x FOR PEER REVIEW b = 0.95 > 1/µ = 0.55, making x (k) converge to zero asymptotically. Finally, 38 Theorem of 44 3 3 3 2(iii) is also verified, since, as it can be observed in Figures 6 and 7, all populations are bounded for all time. 300 250 x (k) 1 x 2 (k) x 3 (k) 200 150 100 50 0 0 5 10 15 20 25 30 35 40 Iteration Figure Figure 6. Evolution 6. Evolution of x1x(k), of x1(k), 2(k) xand 2 (k) x 3(k) xfor and 3 (k)constant for constant parameters parameters in Example in Example 4 when 4 when species species 1 1 dominates dominates the other the other two. two. In case, In this this case, therethere is neither is neither harvesting harvesting nor external nor external consumption. consumption. Example 300 5. Let us consider the multi-species Beverton–Holt Equations (2) and (3) parameterized by the constant values b1 (k) = 0.98, b2 (k ) = 0.96 and b3 (k) = 0.95, di (k) = 0, K1 (k) = 300, K2 (k ) = 200 and, Ke (k) = 200 and µ1 (k) = 2.4, µ2 (k) = 2 and µ3 (k ) = 1.8 with x1 (0) = 100, 250 - x2 (0) = 50 and x3 (0) = 30. The coupling matrix isxgiven 1 (k) by x - (k) 2  -  200 1 x 3 (k) 0.004 0.003 C =  0.0025 1 0.004  150 0.005 0.0035 1 100 The evolution of the three species under this parameterization is displayed in Figure 8. It can be seen that the population of the three species converge to a positive equilibrium 50 point given by x1 = 288.41, x2 = 182.23 and x3 = 175.90, in accordance with Theorem 3, since Condition C1 of Theorem 3(ii) holds. Furthermore, these values are close to the ones 0 calculated 0 5from (31) 10 and 15 given 20 by x 251 =288.46, 30 x352 = 182.58 40 and x3 = 175.42, since the elements Iteration small and satisfy the conditions (27). of the matrix C are sufficiently 300 250 x 1 (k) x 2 (k) 200 x 3 (k) 150 100 0 0 5 10 15 20 25 30 35 40 Iteration Appl. Sci. 2021, 11, 9020 Figure 6. Evolution of x1(k), x2(k) and x3(k) for constant parameters in Example 4 when species 33 1 of 38 dominates the other two. In this case, there is neither harvesting nor external consumption. 300 250 x - (k) 1 x - (k) 2 - 200 x 3 (k) 150 100 50 0 0 5 10 15 20 25 30 35 40 Iteration 300 250 x 1 (k) x 2 (k) 200 x 3 (k) 150 100 50 0 0 5 10 15 20 25 30 35 40 Iteration FigureFigure 7. Evolution 7. Evolution of xx1 2(k), of x1(k), (k) xand x3(k)x3for 2 (k) and (k) constant for constant parameters parameters in Example in Example 4. 4. Species Species 1 1domi- dominates nates the other two and harvesting quota and consumption are included. Example 6. Let us consider the multi-species Beverton–Holt Equations (2) and (3) parameterized by the constant values bi (k) = 1, di (k) = 0, Ki (k) = 200 and µ1 (k) = 2.4, µ2 (k ) = 2 and µ3 (k) = 1.8 with x1 (0) = 100, x2 (0) = 50 and x3 (k) = 30. The coupling matrix is given by   1 0.004 0.005 C =  0.0025 1 0.004  0.005 0.0035 1 The evolution of the three species under this parameterization is displayed in Figure 9. It can be observed, in Figure 9, that there exists a consensus equilibrium point in which the equilibrium points of all species converge to the same unique value. The numerical location of the equilibrium point is given by x = 198.50, which is close to the theoretically calculated one of 198.70 predicted by Proposition 4 in (43). The evolution of the three species under this parameterization is displayed in Figure 8. It can be seen that the population of the three species converge to a positive equilibrium point given by x1 = 288.41, x2 = 182.23 and x3 = 175.90, in accordance with Theorem 3, since Condition C1 of Theorem 3(ii) holds. Furthermore, these values are close to the ones cal- Appl. Sci. 2021, 11, 9020 culated from (31) and given by x1 =288.46, x2 = 182.58 and x3 = 175.42, since the elements of 34 of 38 the matrix C are sufficiently small and satisfy the conditions (27). 300 x -1 (k) 250 x -2 (k) x - (k) 3 200 150 100 50 0 0 10 20 30 40 50 60 70 80 Iteration 300 x 1 (k) 250 x 2 (k) x (k) 3 200 150 Appl. Sci. 2021, 11, x FOR PEER REVIEW 40 of 44 100 50 0 10 20 30 40 50 60 70 80 the equilibrium points of all species Iteration converge to the same unique value. The numerical location of the equilibrium point is given by x = 198.50, which is close to the theoretically FigureFigure 8. Evolution of x1(k),ofx2x(k) 8. Evolution 1 (k), andx2 x(k) 3(k) forxconstant and 3 (k) for constant parametersparameters in Examplein5Example 5 when harvesting when harvesting calculated one of 198.70 predicted by Proposition 4 in (43). quotaquota is included and coupling matrix C has small non-diagonal components. is included and coupling matrix C has small non-diagonal components. Example 200 6. Let us consider the multi-species Beverton–Holt Equations (2) and (3) param- eterized 180 by the constant values bi(k) = 1, di(k) = 0, K i(k) = 200 and μ1 ( k ) = 2.4 , μ 2 ( k ) = 2 - x (k) 1 and μ160 3 ( k ) = 1.8 with x1(0) = 100, x2(0) = 50 and x3(k) = 30. The coupling matrix is given by x (k) - 2 - x (k) 3 140  1 0.004 0.005    120 C =  0.0025 1 0.004   0.005 0.0035 1  100  80 The evolution of the three species under this parameterization is displayed in Figure 9. It can60 be observed, in Figure 9, that there exists a consensus equilibrium point in which 40 20 0 5 10 15 20 25 30 Iteration 200 x 1 (k) x (k) 2 x 3 (k) 150 100 50 0 5 10 15 20 25 30 Iteration FigureFigure 9. Evolution of x1(k),of 9. Evolution x2(k) andx2x(k) x1 (k), 3(k)and for constant x3 (k) for parameters in Examplein6 Example constant parameters when no harvesting 6 when no harvesting quota quota nor external consumption nor external are included consumption so that left-hand are included side values so that left-hand sideand right-hand values side val- side values and right-hand ues are coincident. The coupling matrix C has small non-diagonal entries and the three species con- are coincident. The coupling matrix C has small non-diagonal entries and the three species converge verge to the same equilibrium point. to the same equilibrium point. Example 7. Let us consider the multi-species Beverton–Holt Equations (2) and (3) param- eterized by the constant values bi(k) = 1, di(k) = 0, K1(k) = 300, K2(k) = K3(k) = 200 and μ1 ( k ) = 1.5 + 0.4 k , μ 2 ( k ) = 1.45 + 0.35 k and μ 3 ( k ) = 1.35 + 0.3 k with x1(0) = 100, x2(0) = 75 and x3(k) = 50. The coupling matrix is given by Appl. Sci. 2021, 11, 9020 35 of 38 Example 7. Let us consider the multi-species Beverton–Holt Equations (2) and (3) parameterized by the constant values bi (k) = 1, di (k) = 0, K1 (k) = 300, K2 (k) = K3 (k) = 200 and µ1 (k) = 1.5 + 0.4k , µ2 (k) = 1.45 + 0.35k and µ3 (k) = 1.35 + 0.3k with x1 (0) = 100, x2 (0) = 75 and x3 (k) = 50. The coupling matrix is given by   1 0.4 0.005 C =  0.25 1 0.4  0.05 0.035 1 The evolution of the three species under this parameterization is displayed in Figure 10. Appl. Sci. 2021, 11, x FOR PEER REVIEW Since we are in the condition of applying Theorem 4(i), because 0 = di < bi x44 41 of i and j=k `+1 ∏ j=k` [bi (k `+1 + 1 − j)µi (k `+1 − j)] ≥ ρ > 1, then there is no asymptotic extinction of the populations, as Figure 10 shows. Moreover, all the populations converge to an equi- librium point which is asymptotically stable, as established by Theorem 5. In addition, μ 2 ( k ) = 0.45 + 0.35 k < 1 and μ 3 ( k ) = 0.35 + 0.3 k < 1 . As it is claimed in Theorem 5(ii), the Figure 11 displays the populations of the three species when µ1 (k ) = 0.5 + 0.4k < 1,  n (k) = 0.45 threeµ2species + 0.35k asymptotically, extinguish 0.35 +(10.3 < 1 and µ3 (k) =since ) Asj =it1 cisji xclaimed − μki (<0 )1. i (0) < 1 . in Theorem 4(ii), the three species extinguish asymptotically, since (1 − µi (0))∑nj=1 c ji xi (0) < 1. 2.6 1 (k) (k) 2.4 2 3 (k) 2.2 2 1.8 1.6 1.4 1.2 0 5 10 15 20 25 30 35 40 45 50 Iteration 300 250 x (k) 1 x 2 (k) x (k) 3 200 150 100 50 0 5 10 15 20 25 30 35 40 45 50 Iteration Figure Figure 10. Evolution 10. Evolution of x1x(k), of x1(k), 2(k)xand x3(k)x3for 2 (k) and (k)constant for constant parameters parameters in Example in Example 7 when 7 when nono harvesting harvest- quotanor ing quota norexternal externalconsumption consumptionare areincluded includedandandgrowth growthrates ratesare are larger larger than than unity. unity. 1.6 1 (k) (k) 1.4 2 3 (k) 1.2 1 50 0 5 10 15 20 25 30 35 40 45 50 Iteration Appl. Sci. 2021, 11, 9020 36 of 38 Figure 10. Evolution of x1(k), x2(k) and x3(k) for constant parameters in Example 7 when no harvest- ing quota nor external consumption are included and growth rates are larger than unity. 1.6 1 (k) (k) 1.4 2 3 (k) 1.2 1 0.8 0.6 Appl. Sci. 2021, 11, x FOR PEER REVIEW 0.4 42 of 44 0.2 0 2 4 6 8 10 12 14 Iteration 140 x 1 (k) x (k) 120 2 x 3 (k) 100 80 60 40 20 0 0 5 10 15 Iteration Figure Figure 11. Evolution 11. Evolution of x1x(k), of x1(k), 2(k)xand x3(k)x3for 2 (k) and (k) constant for constant parameters parameters in Example in Example 7 when 7 when nono harvesting harvest- quota nor external consumption are included and growth rates are smaller ing quota nor external consumption are included and growth rates are smaller than unity. than unity. In conclusion, In conclusion, this this section section illustrates, illustrates, through through numerical numerical simulation simulation examples, examples, somesome of the theoretically discussed results. of the theoretically discussed results. 6. Conclusions 6. Conclusions This paper discusses extinction and non-extinction conditions obtained from an This paper discusses extinction and non-extinction conditions obtained from an im- impulsive-type competition Beverton–Holt equation, which, in the most general case, pulsive-type competition Beverton–Holt equation, which, in the most general case, is is modeled under a time-varying parameterization. Some of the obtained results rely modeled under a time-varying parameterization. Some of the obtained results rely on the on the existence of one or more competing species within the studied habitat having a existence of one or more competing species within the studied habitat having a best fit- best fitness, which is interpreted as a dominance of its population stock. Such mentioned ness,extinction/non-extinction which is interpreted as a conditions dominanceare of its population easily testable stock. Such on constraints mentioned extinc- typ- the harvesting, tion/non-extinction conditions are easily testable constraints on the harvesting, ically being hunting/fishing quotas, or, alternatively, on the independent consumption, typically being hunting/fishing which quotas, are the relevant or, alternatively, parameters on the independent in the discontinuities consumption, at sampling which time instants of the are the relevant parameters in the discontinuities at sampling time instants population dynamics. The performed research includes the study of the existence of of the popula- tion extinction dynamics.and The non-extinction performed research includes equilibrium the study points, of the existence the conditions of extinction and of non-negativity and boundedness non-extinctionofequilibrium points, the conditions of non-negativity and boundedness the solutions for any given set of finite non-negative initial conditions, as of the solutions well forgeneric as formal any given set of results finite to related non-negative initialconditions sufficiency-type conditions, forasasymptotic well as formal stability generic results related to sufficiency-type conditions for asymptotic stability to any equilibrium point and, in particular, for extinction. The case of intrinsic to any equi-rates growth librium beingpoint lessand, thaninunity particular, for extinction. and implying Theiscase extinction alsoofdiscussed intrinsic growth compared rates tobeing less the evolution thanofunity and implying the independent extinction isand consumptions alsothediscussed harvesting compared to the evolution quotas sequences. Several of the simulated independent consumptions and the examples are displayed and discussed. harvesting quotas sequences. Several simulated ex- amples are displayed and discussed. Author Contributions: Formal analysis, M.D.l.S. and S.A.-Q.; supervision, S.A.-Q., A.I. and A.J.G.; Software and pro-gramming, A.I.; writing original draft, M.D.l.S., S.A.-Q., A.I, I.G. and A.J.G.; fun- ding acquisition, M.D.l.S., A.J.G. and I.G.; investigation, M.D.l.S., S.A.-Q., A.I, I.G. and A.J.G. All Appl. Sci. 2021, 11, 9020 37 of 38 Author Contributions: Formal analysis, M.D.l.S. and S.A.-Q.; supervision, S.A.-Q., A.I. and A.J.G.; Software and pro-gramming, A.I.; writing original draft, M.D.l.S., S.A.-Q., A.I., I.G. and A.J.G.; funding acquisition, M.D.l.S., A.J.G. and I.G.; investigation, M.D.l.S., S.A.-Q., A.I., I.G. and A.J.G. All authors have read and agreed to the published version of the manuscript. Funding: This research was supported by the Spanish Government through grant RTI2018-094336-B- 100 (MCIU/AEI/FEDER, UE) and by the Basque Government through grant IT1207-19. Institutional Review Board Statement: No applicable. Informed Consent Statement: No applicable. Data Availability Statement: No applicable. Acknowledgments: The authors are grateful to the Spanish Government for its support through grant RTI2018-094336-B-100 (MCIU/AEI/FEDER, UE) and to the Basque Government for its support through grant IT1207-19. The authors are also grateful to the reviewers for their useful comments. Conflicts of Interest: The authors declare that they do not have any conflict of interest. References 1. Barrowman, N.J.; Myers, R.A.; Hilborn, R.; Kehler, D.G.; Field, C.A. The variability among populations of coho salmon in the maximum productive rate and depensation. Ecol. Appl. 2003, 13, 784–793. [CrossRef] 2. Jensen, A.L. Harvest reference points for the Beverton and Holt dynamic pool model. Fish. Res. 2000, 47, 93–96. [CrossRef] 3. Holden, M. Beverton and Holt revisited. Fish. Res. 1995, 24, 3–8. [CrossRef] 4. Baldursson, F.M.; Danielsson, A.; Stefansson, G. On the rational utilization of the Icelandic cod stock. ICES J. Mar. Sci. 1996, 53, 643–658. [CrossRef] 5. Beverton, R.J.H.; Holt, S.J. On the dynamics of exploited fish populations. Fish. Invest. 1957, 19, 1. 6. Elaydi, S.; Sacker, R.J. Global stability of periodic orbits of non-autonomous difference equations and population biology. J. Differ. Equ. 2005, 208, 258–273. [CrossRef] 7. Kocic, V.L. A note on the nonautonomous Beverton-Holt model. J. Differ. Equ. Appl. 2005, 11, 415–422. [CrossRef] 8. Cushing, J.M.; Henson, S.M. A periodically forced Beverton-Holt equation. J. Differ. Equ. Appl. 2002, 8, 1119–1120. [CrossRef] 9. Cushing, J.M.; Henson, S.M. Global dynamics of some periodically forced monotone difference equations. J. Differ. Equ. Appl. 2021, 7, 859–872. [CrossRef] 10. Elaydi, S.; Sacker, R.J. Nonautonomous Beverton-Holt equations and the Cushing-Henson conjectures. J. Differ. Equ. Appl. 2005, 11, 337–346. [CrossRef] 11. Elaydi, S.; Sacker, R.J. Periodic difference equations, population biology and the Cushing-Henson conjectures. Math. Biosci. 2006, 201, 195–207. [CrossRef] 12. Stevic, S. A short proof of the Cushing-Henson conjecture. Discret. Dyn. Nat. Soc. 2006, 2006, 37264. [CrossRef] 13. Berezansky, L.; Braverman, E. On impulsive Beverton-Holt difference equations and their applications. J. Differ. Equ. Appl. 2004, 10, 851–868. [CrossRef] 14. Elabbasy, E.M.; El-Metwally, H.A.; Elsayed, E.M. On the difference equation X(N + 1) = AX(N) − ((BX(N))/(CX(N) − DX(N − 1))). Adv. Differ. Equ. 2006. [CrossRef] 15. de la Sen, M.; Alonso-Quesada, S. Model matching via multirate sampling with fast sampled input guaranteeing the stability of the plant zeros. IEE Proc. Control. Theory Appl. 2007, 1, 210–225. [CrossRef] 16. McCarthy, M.A. The Allee effect, finding mates and theoretical models. Ecol. Model. 1997, 103, 99–102. [CrossRef] 17. Cang, H. Carrying capacity, population equilibrium, and environment’s maximal load. Ecol. Model. 2006, 192, 317–320. 18. de la Sen, M. The environment carrying capacity is not independent of the intrinsic growth rate for subcritical spawning stock biomass in the Beverton-Holt equation. Ecol. Model. 2007, 204, 271–273. [CrossRef] 19. Meng, X.; Xu, W.; Chen, L. Profitless delays for a nonautonomous Lotka-Volterra predator–prey almost periodic system with dispersion. Appl. Math. Comput. 2007, 188, 365–378. [CrossRef] 20. White, L.; White, F.; Luo, Y.; Xu, T. Estimation of parameters in carbon sequestration models from net ecosystem exchange data. Appl. Math. Comput. 2006, 181, 864–879. [CrossRef] 21. Sarkar, R.R.; Mukhopadhyay, B.; Bhattacharyya, R.; Banerjee, S. Time lags can control algal bloom in two harmful phytoplankton– zooplankton system. Appl. Math. Comput. 2007, 186, 445–459. [CrossRef] 22. Chen, F.; Shi, C. Global attractivity in an almost periodic multi-species nonlinear ecological model. Appl. Math. Comput. 2006, 180, 376–392. [CrossRef] 23. Villacampa-Esteve, Y.; Castro-Lopez, M.A. A text theory of ecological systems. Cybern. Syst. 1999, 30, 587–607. [CrossRef] 24. Usó-Domènech, J.L.; Villacampa-Esteve, Y.; Mateu-Mahiques, J.; Sastre-Vazquez, P. Uncertainty and complementarity principles in ecological models. Cybern. Syst. 2000, 31, 137–159. 25. de la Sen, M.; Luo, N.S. A note on the stability of linear time-delay systems with impulsive inputs. IEEE Trans. Circuits Syst. I.Fundam. Theory Appl. 2003, 50, 149–152. [CrossRef] Appl. Sci. 2021, 11, 9020 38 of 38 26. de la Sen, M.; Alonso-Quesada, S. On the impulsive Beverton-Holt equation and extinction conditions in population dynamics. Nat. Resour. Modeling 2011, 24, 514–534. [CrossRef] 27. Shalsh, O.K.; Al-Nassir, S. Dynamics and optimal harvesting strategy for biological models with Beverton-Holt growth. Iraqi J. Sci. 2020, 223–232. [CrossRef] 28. de la Sen, M.; Alonso-Quesada, S. A control theory point of view on Beverton-Holt equation in population dynamics and some of its generalizations. Appl. Math. Comput. 2008, 199, 464–481. [CrossRef] 29. de la Sen, M.; Alonso-Quesada, S. Model-matching based control of the Beverton-Holt equation in Ecology. Discret. Dyn. Nat. Soc. 2008, 2008, 793512. [CrossRef] 30. Bertrand, E.; Kulenovic, M.R.S. Global dynamics of generalized second-order Beverton-Holt equations of linear and quadratic type. J. Comput. Anal. Appl. 2021, 29, 185–202. 31. Khyat, T.; Kulenovic, M.R.S. Global dynamics of delayed sigmoid Beverton-Holt equation. Discret. Dyn. Nat. Soc. 2020, 2020, 1364282. [CrossRef] 32. Al-Ghassani, A.S.; Alsharawi, Z. The effect of maps permutation on the global attractor of a periodic Beverton-Holt model. Appl. Math. Comput. 2020, 370, 124905. [CrossRef] 33. de la Sen, M. The generalized Beverton-Holt equation and the control of populations. Appl. Math. Comput. 2008, 32, 2312–2328. [CrossRef] 34. Kang, Y. Dynamics of a generalized Ricker-Beverton-Holt competition model subject to Allee effects. J. Differ. Equ. Appl. 2016, 22, 687–723. [CrossRef] 35. de la Sen, M.; Ibeas, A.; Garrido, A.J. Stage-dependent structured discrete-time models for mosquito population evolution with survability: Solution properties, equilibrium points, oscillations, and populatio feedback controls. Mathematics 2019, 7, 1181. [CrossRef]