# Principle of Limiting Factors-Driven Piecewise Population Growth Model I: Qualitative Exploration and Study Cases on Continuous-Time Dynamics

Authors: Héctor A. Echavarria-Heras; Cecilia Leal-Ramírez; Guillermo Gómez; Elia Montiel-Arzate
Journal: Complexity
 (2021)

---

## Abstract

We examine the comportment of the global trajectory of a piecewisely conceived single species population growth model. Formulation relies on what we develop as the principle of limiting factors for population growth, adapted from the law of the minimum of Liebig and the law of the tolerance of Shelford. The ensuing paradigm sets natality and mortality rates to express through extreme values of population growth determining factor. Dynamics through time occur over different growth phases. Transition points are interpreted as thresholds of viability, starvation, and intraspecific competition. In this delivery, we focus on the qualitative study of the global trajectory expressed on continuous time and on exploring the feasibility of analytical results against data on populations growing under experimental or natural conditions. All study cases sustained fittings of high reproducibility both at empirical and interpretative slants. Possible phase configurations include regimes with multiple stable equilibria, sigmoidal growth, extinction, or stationarity. Here, we also outline that the associating discrete-time piecewise model composes the logistic map applied over a particular region of the phase configuration. Preliminary exploratory analysis suggests that the logistic map’s chaos onset could surpass once the orbit enters a contiguous phase region.


---

## Body

## 1. Introduction

As a rule, mathematical models that attempt to describe the dynamics of two or more populations, subject to a specific type of interaction in an ecosystem, are formulated based on too rigid hypotheses, resulting in poor concordance of the model’s predictions with the natural process. The high degree of complexity of biological phenomena and their significant spatial-temporal variations do not assimilate by single-rule models, which are incapable of presenting the functional diversity required by levels of reliable prediction. However, when trying to include more biological information in the model hypotheses, it is often possible to fall into the opposite situation: creating models or techniques so complicated that they are also inoperative from an analytical standpoint.From the traditional point of view, in building a model, it is necessary to include the most representative processes of the system that we are trying to describe to maintain simplicity as much as possible without losing relevant information. The problem reduces to a quest for those key processes that govern the interaction dynamics we are interested in studying. For example, in demographic ecology, the dynamics of a population depend on the correlation between opposing influences such as birth generating and growth inhibitory processes. Each of these effects depends on a series of factors inherent to population growth. For example, births depend on age composition, sex proportion, fertility, available food, and so on. Correspondingly, the growth inhibition process may be set by population density, disease, abiotic factors such as temperature, humidity, pollution, and others. As the number of determining factors increases, the model’s complexity affects tractability and interpretability. An alternative approach is considering a horizontal integration of complexity through a piecewise modeling strategy. While embracing this approach to gain interpretative strength, it would be desirable to imbue the construct of a mechanistic profile. This work explains how we build a piecewise population model (PPM) by following a systematic way, that is, creating a piecewise description of population growth by relying upon a logical deductive approach. An ad hoc logical system that could sustain such an enterprise is what we comprehend as a principle of limiting factors (PLF) for population growth. This paradigm entailed the derivation of a PPM as a collection of submodels called growth phases that continuously compose the global trend of population size. Each growth phase describes the dynamics over a specific time interval. The naming composite further refers to the principle of limiting factors-driven piecewise population model, which also represents utilizing the PLF-PPM acronym.As we elucidate in what follows, the PLF-PPM can be arranged by following ideas in [1–3]. Mainly, the construct presented here reviews the approach in [3] to consider a modification taking into account the abatement effect that the population induces on an external resource necessary to guarantee its growth and permanence. The present PLF-PPM formulation also includes a specific scaling or weighting of population size to model the increase in mortality promoted by low population densities [4]. The resulting device allows meaningful ecological interpretation of consequential growth phases. Mainly for keeping this paper’s extension manageable, we focus on the continuous-time form of the PLF-PPM. Nevertheless, we judged that it is pertinent to advance an outline of the next discrete-time version, whose detailed exploration will be addressed in the second part of this work. We include several study cases that show the empirical and interpretative adequacy of the present paradigm. The Appendix presents a qualitative study of associating continuous-time global trajectory.


## 2. Materials and Methods

### 2.1. The General Piecewise Population Growth Model Setup

Throughout this paper, symbolR stands for the set of real numbers and R+ denotes the subset of R whose elements are positive. Correspondingly, the size of a population at a positive time t represents employing xt and formally stands as a function having domain R+ and range XR, a proper subset of R+, that is, xt:R+⟶XR⊆R+. For the case of an isolated population, a customary assumption establishes the existence of a function Gx:XR⟶R such that(1)x˙t=Gxt⋅xt,where x˙t stands for the rate of change of the number xt of individuals in the population at time t. The function Gx is interpreted as an average estimate of the influence of all the processes that govern population growth. More specifically, this function embodies the average growth rate for each individual in the population. In the following, we will refer to Gt as the per capita growth rate or the intrinsic growth rate of xt.Determining the functionGx in a closed form is usually impossible or very complicated. However, from an empirical perspective, it is possible to suppose that the function Gx admits expansion as a power series, with which we would obtain(2)Gx=∑i=0∞aixi.Then, for example, neglecting theanxn terms with n≥2, we could consider the approximation(3)Gx=r−kx.The above characterization of the intrinsic growth rateGx establishes the form of the differential equation (1) that leads to the identification of the so-called logistic curve, commonly adopted as a model of population growth in a limited environment. Although it is feasible to give a biological interpretation to equation (3) model, it shows the fragility inherent in all simplified models, which mainly relates to the lack of suitable complexity. Moreover, the variation of population size could compose an outstanding array of patterns. Then, maintaining a single association rule for Gx, through the full extension of XR,no matter how complex such a rule conceives, it could render insufficient to imbue a consistent reproducibility of observed population values. A traditional approach to adapt Gx is by considering a polynomial approximation of degree n derived from the infinite series of equation (2). Nevertheless, this approach bears a disadvantage since at gaining reproducibility strength, we could lose interpretability.Before getting into matters on suitable complexity embedding ofGx, we introduce some notation conventions. Given a functionSt:A⟶B and a subset C⊆A, the restriction of St to the set, C denotes St|C:C⟶B, it is formally defined by (4)St|C=Stfor eachtin C.A second method of adapting the complexity attending to an optimal reproducibility criterion is conceivingGx as a composite of m operating modes, each one characterized by a continuous and differentiable function Gix:ℜix⟶R being ℜix a region defined through(5)ℜix=x|pi−1<x≤pi,for i=1,2,…,m. The components of the vector P=p0,…,pm represent breakpoints for transition among the operating modes Gix. We also assume that the collection of regions ℜix provides a covering for the range of xt, that is, XR⊆∪1mℜix. We could formally set(6)Gix=Gx|ℜix.So, we can representGxexpanded by the Gix submodels for i=1,2,…,m, namely,(7)Gx=∑i=1mχℜixGix,with χℜix representing the characteristic function of ℜix, namely,(8)χℜix=1,ifxt∈ℜix,0,ifxt∈ℜixc,1.Equation (7) intends to provide a piecewise description of the intrinsic population growth rate in the autonomous differential equation (1). Discretizing the resulting differential equation, we can identify the parameters defining the Gix submodels describing the phases of population growth, as well as the threshold parameters p0,p1,…,pm that represent the transition points between the growth phases Gix; given assumed continuity properties of Gx, the phase changes must be set as continuous transitions between the submodels, that is, we have(9)Gipi=limx⟶pi−Gi+1x,fori=1,2,…,m−1.But, in the general settings, the transition between phases at breakpointsp0,p1,…,pm can occur in a discontinuous way.In the settings of the model of equation (1), the function fx=Gxx is interpreted as the natural growth rate of the population. Then, defining fix=Gixx for i=1,2,…,m and keeping the notation convention entailing the associated domain regions, according to equation (6), we set(10)fix=fx|ℜix.This way, we can propose a piecewise expanded form of equation (1), namely,(11)x˙t=∑i=1mχℜixfix,forx0=x0,where ℜix and χℜix are defined by equations (5) and (8) one to one.The non-autonomous form of the dynamical system of equation (11) addresses similarly by replacing fx by a continuous and differentiable function fx,t:XR×R⟶R and interprets xt and t as a forced form of the natural growth rate of the population. The extension of equation (11) to the non-autonomous case fx,t follows by setting submodel fix,t:ℜix,t⟶R being ℜix,t=ℜix×R for i=1,2,…,m. This way, a piecewisely expanded non-autonomous form of equation (11) could be expressed by(12)x˙t=∑1=1mχℜix,tfix,t,withx0=x0,where χℜix,t symbolizes the characteristic function of the rectangle ℜix,t. Nevertheless, in this paper, we focus on the autonomous case.To provide a piecewise representation of the solutionxt to the dynamical system of equation (11), we conceive trajectory sectors xit for i=1,2,…,m defined by(13)xit=xt|Tit.ForTit⊆R+, the interval Tit=x−1ℜix. Since xt varies continuously, we must set the consistency condition(14)xiti=limt⟶ti−xi+1t,for ti=x−1pi being pi the breakpoint for the transition from fix to fi+1x. Then, we arrange(15)xt=∑i=1mχTitxit,forx0=x0,with χTit being defined by(16)χTit=1,ift∈Tit,0,ift∈Titc.In any event, usually assembling equation (11) relies on an empirical description that abides by the highest reproducibility criterion. Regularly, such an approach disregards any phenomenological explanation of population growth. In this work, we explain a systematic way to obtain a collection of submodels that mimic the global trend of population size xt, where each one of them describes the dynamics over a specific time interval, that is, we build a piecewise description of both x˙t and xt as given by equations (11) and (13) one to one by relying on a logical deductive approach. This can be obtained by following [1–3]. The corresponding derivation is based on the use of the ecological principle of limiting factors. The construct presented here corresponds to a modification of the slant in [3]. To consider the abatement effect that the population induces on an external resource necessary to guarantee its growth and permanence, the present formulation also includes a specific scaling or weighting of population size to model the increase in mortality promoted by low population densities [4].


### 2.2. The Principle of Limiting Factors

German Physiologist Justus Von Liebig [5], studying the growth of certain plants, realized that to guarantee their development, there had to be a set of essential nutrients. Some had to be abundant, and others were required only in small quantities. A significant discovery of Liebig was that the absence of some nutrient could not be replaced with any other that appeared in abundance. Moreover, a medium that contained all the nutrients in abundance except one of them, which appeared in insufficient quantity, would allow plant growth until the nutrient was utterly depleted. The absence of some element limited growth. In addition, when growth occurred, the latter turned out to be controlled by the nutrient appearing in lesser proportion. Liebig called this regularity the law of the minimum. Later, in [6], the North American Ecologist V. E. Shelford extended the domain of Liebig’s law to what is now known as the law of tolerance. Shelford pointed out that when there is an excess of a specific element, this can be as limiting as its deficiency. It follows that all the processes that determine the dynamics of a population will occur in intensities governed by the minima or maxima of the factors that nourish these processes. We refer to this notion as the principle of limiting factors (PLF) for simplicity and following the tradition. It intends to interpret the law that limits a process through the maximum and minimum of factors acting in the analyzed system.


### 2.3. The Principle of Limiting Factors-Driven Piecewise Population Growth Model

Considering the validity of the principle of limiting factors, let us construct a mathematical model to interpret the performance of a population whose size at timetis denoted by xt. We suppose that departing from a value x0=x0, the population growth depends on an external resource whose availability at time t is modeled by a continuous and differentiable function Rt with rank in R+. We will also assume that the natural growth rate of the population size x˙t as given by equation (2) expresses through a function fx=Gxx that determines on a first instance by the intensities of two opposing processes, one of birth Nx yielding offspring and another Mx determining the death of the individuals, that is,(17)fx=βNx−μMx,where β and μ are positive scaling parameters [2, 3], Nx stands for the strength of the process leading to increasing population size, that is, the rate at which the number of newborn in the population is generated, and Mx represents the intensity of the growth-inhibiting process, that is, the mortality rate. We assume that Nx and Mx both express through continuous and differentiable functions.Then, according to the principle of limiting factors, we assume thatNx depends essentially on two factors: one being the size of the population xt and the second one taken as the amount of feeding energy available; this symbolized employing Rx. We suppose that before the population began exhausting the resource Rx, this is kept at a constant value F0. Additionally, we think through that the total weight of the population expresses in the same units of energy as Rx. Let us also undertake that each individual in the xt pool maintains an energy consumption equal to a proportion k of its weight. Then, we can assume in the first instance that Rx represents using the linear model(18)Rx=F0−kxt.Then, resource availabilityRxo will plentifully suit feeding requirements of the initial population xo whenever the condition Rxo≫kxo satisfies or equivalently according to equation (18) whenever F0≫2kx0 holds. We will then have that if the amount of energy consumed by the population at time t, i.e., kxt, is less than Rx, there will be enough food for all individuals, that is, the availability of energy will guarantee that the population can participate in the reproduction process unrestrictedly, so we can establish the assumption that the birth rate will be proportional to the size of the population, that is, we can set the equality Nx=xt. In the case in which kxt exceeds Rx, then only as many individuals of the population energetically equivalent to Rx will be able to feed. Therefore, the reproductive potential of the population will be limited. That is, according to Liebig’s principle, energy availability will be a limiting factor for population growth. Then, Nx=Rx. Thus, we can define Nx by(19)Nx=mintkxt,Rx.Solving the inequalitykxt<F0−kxt, we can express equation (19) in the equivalent form(20)Nx=mintxt,E,where the constant E=F0/2k is interpreted as a starvation threshold since if xt grows maintaining levels higher than E, the number of offspring will turn to be controlled by the external factor F0 and the rate of individual energy consumption k that could eventually lead to the decline of the population.Concerning the intensity of the natural death processMx, we will consider this, considering Shelford’s law of tolerance based on population density as the critical factor. According to Allee’s principle, reduced population levels will restrict the number of interactions between individuals, making mating difficult, implying decreasing recruitment, favoring even a lower number of individuals, and eventually leading to the total depletion of the population. In this circumstance, we will say that the strength of the mortality process will be proportional to the size of the population, that is, Mx=axt, where a is a positive constant. On the other hand, at high population sizes, density upsurges, the struggle for resources intensifies, and so does the transmission of diseases, thereby increasing mortality. In this case, we take the intensity of the mortality process as being proportional to the number of encounters between the individuals, that is, Mx=bx2t where b∈R+. Then, according to Shelford’s law of tolerance, we can express Mx in the form(21)Mx=maxtaxt,bx2t.In summary, the PLF and equations (11) and (17) imply the piecewise continuous form(22)x˙t=∑i=1mχℜixfix,forx0=x0,where fix and ℜixi=1,2,3,4 are given by(23)f1x=βE−μbxt2ℜ1x=xt|xt>E∧xt>abf2x=βxt−μbxt2ℜ2x=xt|xt<E∧xt>abf3x=β−μaxtℜ3x=xt|xt<E∧xt<abf4x=βE−μaxtℜ4x=xt|xt>E∧xt<ab,with χℜix representing the characteristic function of ℜix for i=1,2,3,4 (cf. equation (8)). In what follows, the principle of limiting factors-driven piecewise population growth model of equations (22) and (23) will be referred to as PLF-PPM for short.Sincext stands for the solution to equation (22), then according to equation (15), we have that xit=xt|Tit for i=1,2,3,4, stands for the ith sector of the PLF driven trajectory xt having range in ℜix, and correspondingly, this sets Tit=x−1ℜix. Then, xit solves the differential equation x˙it=fixt with xit satisfying xi0=xi0∈ℜix and the restriction t∈Tit. Therefore, we have(24)x1t=λ1C−e−2tβμbEC+e−2tβμbExt∈ℜ1x,t∈T1tx10=λ1C−1C+1λ1=βEμbx2t=Cλ2λ2e−βt+Cxt∈ℜ2x,t∈T2tx20=βCβ+μbCλ2=βμbx3t=Ceβ−μatxt∈ℜ3x,t∈T3tx30=Cx4t=Ce−μat+λ41−e−μatxt∈ℜ4x,t∈T4tx40=Cλ4=βEμa,with λi being the equilibrium solutions as determined by the natural growth rates fixt (cf. equation (24)). Additionally, the submodels xit composing the solution xt as given by equation (15) determine by the order relationship that E and a/b satisfy. Accordingly, the global trajectory xtcould fit to three possible topologies, one type M1associating to the statement M1:E>a/b, another form M2 consistent to M2:E>a/b; and a third form M3 associating to M3:E>a/b the phase composite P1t,xt for M1 models will include regions ℜ1x,ℜ2x, and ℜ3x, that is, P1t,xt=∪i≠4ℜix∪​0,∞. Correspondingly, for models M2, phase arrangement becomes P2t,xt=∪i≠2ℜix∪​0,∞. In turn, phase portrait for models M2 and the regions ℜ2x and ℜ4x is mutually exclusive. For models M3, the phase combination turns to P3t,xt=∪i≠2,4ℜix∪​0,∞. The Appendix presents a qualitative study of trajectories resulting from equation (22).


## 2.1. The General Piecewise Population Growth Model Setup

Throughout this paper, symbolR stands for the set of real numbers and R+ denotes the subset of R whose elements are positive. Correspondingly, the size of a population at a positive time t represents employing xt and formally stands as a function having domain R+ and range XR, a proper subset of R+, that is, xt:R+⟶XR⊆R+. For the case of an isolated population, a customary assumption establishes the existence of a function Gx:XR⟶R such that(1)x˙t=Gxt⋅xt,where x˙t stands for the rate of change of the number xt of individuals in the population at time t. The function Gx is interpreted as an average estimate of the influence of all the processes that govern population growth. More specifically, this function embodies the average growth rate for each individual in the population. In the following, we will refer to Gt as the per capita growth rate or the intrinsic growth rate of xt.Determining the functionGx in a closed form is usually impossible or very complicated. However, from an empirical perspective, it is possible to suppose that the function Gx admits expansion as a power series, with which we would obtain(2)Gx=∑i=0∞aixi.Then, for example, neglecting theanxn terms with n≥2, we could consider the approximation(3)Gx=r−kx.The above characterization of the intrinsic growth rateGx establishes the form of the differential equation (1) that leads to the identification of the so-called logistic curve, commonly adopted as a model of population growth in a limited environment. Although it is feasible to give a biological interpretation to equation (3) model, it shows the fragility inherent in all simplified models, which mainly relates to the lack of suitable complexity. Moreover, the variation of population size could compose an outstanding array of patterns. Then, maintaining a single association rule for Gx, through the full extension of XR,no matter how complex such a rule conceives, it could render insufficient to imbue a consistent reproducibility of observed population values. A traditional approach to adapt Gx is by considering a polynomial approximation of degree n derived from the infinite series of equation (2). Nevertheless, this approach bears a disadvantage since at gaining reproducibility strength, we could lose interpretability.Before getting into matters on suitable complexity embedding ofGx, we introduce some notation conventions. Given a functionSt:A⟶B and a subset C⊆A, the restriction of St to the set, C denotes St|C:C⟶B, it is formally defined by (4)St|C=Stfor eachtin C.A second method of adapting the complexity attending to an optimal reproducibility criterion is conceivingGx as a composite of m operating modes, each one characterized by a continuous and differentiable function Gix:ℜix⟶R being ℜix a region defined through(5)ℜix=x|pi−1<x≤pi,for i=1,2,…,m. The components of the vector P=p0,…,pm represent breakpoints for transition among the operating modes Gix. We also assume that the collection of regions ℜix provides a covering for the range of xt, that is, XR⊆∪1mℜix. We could formally set(6)Gix=Gx|ℜix.So, we can representGxexpanded by the Gix submodels for i=1,2,…,m, namely,(7)Gx=∑i=1mχℜixGix,with χℜix representing the characteristic function of ℜix, namely,(8)χℜix=1,ifxt∈ℜix,0,ifxt∈ℜixc,1.Equation (7) intends to provide a piecewise description of the intrinsic population growth rate in the autonomous differential equation (1). Discretizing the resulting differential equation, we can identify the parameters defining the Gix submodels describing the phases of population growth, as well as the threshold parameters p0,p1,…,pm that represent the transition points between the growth phases Gix; given assumed continuity properties of Gx, the phase changes must be set as continuous transitions between the submodels, that is, we have(9)Gipi=limx⟶pi−Gi+1x,fori=1,2,…,m−1.But, in the general settings, the transition between phases at breakpointsp0,p1,…,pm can occur in a discontinuous way.In the settings of the model of equation (1), the function fx=Gxx is interpreted as the natural growth rate of the population. Then, defining fix=Gixx for i=1,2,…,m and keeping the notation convention entailing the associated domain regions, according to equation (6), we set(10)fix=fx|ℜix.This way, we can propose a piecewise expanded form of equation (1), namely,(11)x˙t=∑i=1mχℜixfix,forx0=x0,where ℜix and χℜix are defined by equations (5) and (8) one to one.The non-autonomous form of the dynamical system of equation (11) addresses similarly by replacing fx by a continuous and differentiable function fx,t:XR×R⟶R and interprets xt and t as a forced form of the natural growth rate of the population. The extension of equation (11) to the non-autonomous case fx,t follows by setting submodel fix,t:ℜix,t⟶R being ℜix,t=ℜix×R for i=1,2,…,m. This way, a piecewisely expanded non-autonomous form of equation (11) could be expressed by(12)x˙t=∑1=1mχℜix,tfix,t,withx0=x0,where χℜix,t symbolizes the characteristic function of the rectangle ℜix,t. Nevertheless, in this paper, we focus on the autonomous case.To provide a piecewise representation of the solutionxt to the dynamical system of equation (11), we conceive trajectory sectors xit for i=1,2,…,m defined by(13)xit=xt|Tit.ForTit⊆R+, the interval Tit=x−1ℜix. Since xt varies continuously, we must set the consistency condition(14)xiti=limt⟶ti−xi+1t,for ti=x−1pi being pi the breakpoint for the transition from fix to fi+1x. Then, we arrange(15)xt=∑i=1mχTitxit,forx0=x0,with χTit being defined by(16)χTit=1,ift∈Tit,0,ift∈Titc.In any event, usually assembling equation (11) relies on an empirical description that abides by the highest reproducibility criterion. Regularly, such an approach disregards any phenomenological explanation of population growth. In this work, we explain a systematic way to obtain a collection of submodels that mimic the global trend of population size xt, where each one of them describes the dynamics over a specific time interval, that is, we build a piecewise description of both x˙t and xt as given by equations (11) and (13) one to one by relying on a logical deductive approach. This can be obtained by following [1–3]. The corresponding derivation is based on the use of the ecological principle of limiting factors. The construct presented here corresponds to a modification of the slant in [3]. To consider the abatement effect that the population induces on an external resource necessary to guarantee its growth and permanence, the present formulation also includes a specific scaling or weighting of population size to model the increase in mortality promoted by low population densities [4].


## 2.2. The Principle of Limiting Factors

German Physiologist Justus Von Liebig [5], studying the growth of certain plants, realized that to guarantee their development, there had to be a set of essential nutrients. Some had to be abundant, and others were required only in small quantities. A significant discovery of Liebig was that the absence of some nutrient could not be replaced with any other that appeared in abundance. Moreover, a medium that contained all the nutrients in abundance except one of them, which appeared in insufficient quantity, would allow plant growth until the nutrient was utterly depleted. The absence of some element limited growth. In addition, when growth occurred, the latter turned out to be controlled by the nutrient appearing in lesser proportion. Liebig called this regularity the law of the minimum. Later, in [6], the North American Ecologist V. E. Shelford extended the domain of Liebig’s law to what is now known as the law of tolerance. Shelford pointed out that when there is an excess of a specific element, this can be as limiting as its deficiency. It follows that all the processes that determine the dynamics of a population will occur in intensities governed by the minima or maxima of the factors that nourish these processes. We refer to this notion as the principle of limiting factors (PLF) for simplicity and following the tradition. It intends to interpret the law that limits a process through the maximum and minimum of factors acting in the analyzed system.


## 2.3. The Principle of Limiting Factors-Driven Piecewise Population Growth Model

Considering the validity of the principle of limiting factors, let us construct a mathematical model to interpret the performance of a population whose size at timetis denoted by xt. We suppose that departing from a value x0=x0, the population growth depends on an external resource whose availability at time t is modeled by a continuous and differentiable function Rt with rank in R+. We will also assume that the natural growth rate of the population size x˙t as given by equation (2) expresses through a function fx=Gxx that determines on a first instance by the intensities of two opposing processes, one of birth Nx yielding offspring and another Mx determining the death of the individuals, that is,(17)fx=βNx−μMx,where β and μ are positive scaling parameters [2, 3], Nx stands for the strength of the process leading to increasing population size, that is, the rate at which the number of newborn in the population is generated, and Mx represents the intensity of the growth-inhibiting process, that is, the mortality rate. We assume that Nx and Mx both express through continuous and differentiable functions.Then, according to the principle of limiting factors, we assume thatNx depends essentially on two factors: one being the size of the population xt and the second one taken as the amount of feeding energy available; this symbolized employing Rx. We suppose that before the population began exhausting the resource Rx, this is kept at a constant value F0. Additionally, we think through that the total weight of the population expresses in the same units of energy as Rx. Let us also undertake that each individual in the xt pool maintains an energy consumption equal to a proportion k of its weight. Then, we can assume in the first instance that Rx represents using the linear model(18)Rx=F0−kxt.Then, resource availabilityRxo will plentifully suit feeding requirements of the initial population xo whenever the condition Rxo≫kxo satisfies or equivalently according to equation (18) whenever F0≫2kx0 holds. We will then have that if the amount of energy consumed by the population at time t, i.e., kxt, is less than Rx, there will be enough food for all individuals, that is, the availability of energy will guarantee that the population can participate in the reproduction process unrestrictedly, so we can establish the assumption that the birth rate will be proportional to the size of the population, that is, we can set the equality Nx=xt. In the case in which kxt exceeds Rx, then only as many individuals of the population energetically equivalent to Rx will be able to feed. Therefore, the reproductive potential of the population will be limited. That is, according to Liebig’s principle, energy availability will be a limiting factor for population growth. Then, Nx=Rx. Thus, we can define Nx by(19)Nx=mintkxt,Rx.Solving the inequalitykxt<F0−kxt, we can express equation (19) in the equivalent form(20)Nx=mintxt,E,where the constant E=F0/2k is interpreted as a starvation threshold since if xt grows maintaining levels higher than E, the number of offspring will turn to be controlled by the external factor F0 and the rate of individual energy consumption k that could eventually lead to the decline of the population.Concerning the intensity of the natural death processMx, we will consider this, considering Shelford’s law of tolerance based on population density as the critical factor. According to Allee’s principle, reduced population levels will restrict the number of interactions between individuals, making mating difficult, implying decreasing recruitment, favoring even a lower number of individuals, and eventually leading to the total depletion of the population. In this circumstance, we will say that the strength of the mortality process will be proportional to the size of the population, that is, Mx=axt, where a is a positive constant. On the other hand, at high population sizes, density upsurges, the struggle for resources intensifies, and so does the transmission of diseases, thereby increasing mortality. In this case, we take the intensity of the mortality process as being proportional to the number of encounters between the individuals, that is, Mx=bx2t where b∈R+. Then, according to Shelford’s law of tolerance, we can express Mx in the form(21)Mx=maxtaxt,bx2t.In summary, the PLF and equations (11) and (17) imply the piecewise continuous form(22)x˙t=∑i=1mχℜixfix,forx0=x0,where fix and ℜixi=1,2,3,4 are given by(23)f1x=βE−μbxt2ℜ1x=xt|xt>E∧xt>abf2x=βxt−μbxt2ℜ2x=xt|xt<E∧xt>abf3x=β−μaxtℜ3x=xt|xt<E∧xt<abf4x=βE−μaxtℜ4x=xt|xt>E∧xt<ab,with χℜix representing the characteristic function of ℜix for i=1,2,3,4 (cf. equation (8)). In what follows, the principle of limiting factors-driven piecewise population growth model of equations (22) and (23) will be referred to as PLF-PPM for short.Sincext stands for the solution to equation (22), then according to equation (15), we have that xit=xt|Tit for i=1,2,3,4, stands for the ith sector of the PLF driven trajectory xt having range in ℜix, and correspondingly, this sets Tit=x−1ℜix. Then, xit solves the differential equation x˙it=fixt with xit satisfying xi0=xi0∈ℜix and the restriction t∈Tit. Therefore, we have(24)x1t=λ1C−e−2tβμbEC+e−2tβμbExt∈ℜ1x,t∈T1tx10=λ1C−1C+1λ1=βEμbx2t=Cλ2λ2e−βt+Cxt∈ℜ2x,t∈T2tx20=βCβ+μbCλ2=βμbx3t=Ceβ−μatxt∈ℜ3x,t∈T3tx30=Cx4t=Ce−μat+λ41−e−μatxt∈ℜ4x,t∈T4tx40=Cλ4=βEμa,with λi being the equilibrium solutions as determined by the natural growth rates fixt (cf. equation (24)). Additionally, the submodels xit composing the solution xt as given by equation (15) determine by the order relationship that E and a/b satisfy. Accordingly, the global trajectory xtcould fit to three possible topologies, one type M1associating to the statement M1:E>a/b, another form M2 consistent to M2:E>a/b; and a third form M3 associating to M3:E>a/b the phase composite P1t,xt for M1 models will include regions ℜ1x,ℜ2x, and ℜ3x, that is, P1t,xt=∪i≠4ℜix∪​0,∞. Correspondingly, for models M2, phase arrangement becomes P2t,xt=∪i≠2ℜix∪​0,∞. In turn, phase portrait for models M2 and the regions ℜ2x and ℜ4x is mutually exclusive. For models M3, the phase combination turns to P3t,xt=∪i≠2,4ℜix∪​0,∞. The Appendix presents a qualitative study of trajectories resulting from equation (22).


## 3. Results

### 3.1. Analytic Exploration

To study the varied forms of the global trajectoryxt derived from the PLF-PM, we begin by setting the array of possible phase configurations associated with the basic parameters a,b,E,β, and μ (cf. equations (22) and (23)). We depart from the ordering relationship for the external energy E and the ratio of mortality from population density proportions, a/b, that yields a first model type classification Mi,i=1,2,3. We have to consider models M1 associated to the inequality E>a/b, where the range of xt composes the regions: ℜ1x=a/b<E<x, ℜ2x=a/b<x<E, and ℜ3x=x<a/b<E (cf. equation (23)). For model type M2 where the relationship E>a/b holds, the range of xt will be divided into the regions ℜ1x,ℜ3x, and ℜ4x=x<a/b<E (cf. equation (23)). For model type M3 linking to the ordering E=a/b, the range of xt will be divided into the regions ℜ1x and ℜ3x. We can be aware that regions ℜ2x and ℜ4x are mutually excluding. Besides, it also requires classifying possible phase arrangements according to a Birth to Mortality Scaled Ordering BMSi,i=1,2,3, namely, BMS1: when β<μa, BMS2: if β>μa, and BMS3: whenever β=μa. Additionally, each model type Mi with i=1,2 determines possible placements of the initial condition, and each one associated to composing phase regions ℜix with i=1,2,3,4 determines an initial condition ordering ICikiwithi=1,2andki=1,2,…,5, namely, IC11: for (E<xo), IC12: whenever a/b<xo<E, IC13: when xo<a/b, IC14: if xo=E, and IC15: whenever xo=a/b. Correspondingly, IC21: for xo<E, IC22: whenever E<xo<a/b, IC23: when xo>a/b,IC24: if xo=E, and IC25: whenever xo>a/b. The possible phase arrangements typify according to a three-dimensional conjunction operator ℋi,j,ki=Mi∧BMSj∧ICiki for iand ki as specified. Table 1 summarizes phase arrangements for model type M1, and Table 2 summarizes those corresponding to model type M2. In addition to the ℋi,j,ki conjunctions, we must consider the positioning of the equilibrium solutions inside their associated regions. For this reason, in the presentation of results, we could use the Oλi,j indicator for i=1,2,4 designating equilibrium λi and j=1,2,4 describing region ℜix. Then, for instance, Oλ1,1 identifies the placement of equilibrium λ1 inside region ℜ1x; correspondingly, Oλ2,3 labels positioning of λ2 inside region ℜ3x and so on. Yet given Oλi,j, this relates to distinguishing how the initial condition is positioned relative to the equilibrium solution λi. Resulting orderings generically symbolize through ICλi,j for i=1,2,4 associated to the equilibrium solution λi and j=1,2 labeling the ordering of x0 relative to λi. This way, ICλi,1 stands for x0<λi and ICλi,2 stands for x0>λi. Besides described symbols, we could refer to direct inequalities that a given equilibrium satisfies relative to thresholds defining regions.Table 1
Basic phase arrangementsℋ1,j,ki determined by the PLF-PM and phase arrangements for model type M1 associated to the E>a/b ordering.
ℋi,j,kiM1BMSiBMSiICikiICikiℋ1,1,1E>a/bBMS1β<μaIC11E<xoℋ1,1,2E>a/bBMS1β<μaIC12a/b<xo<Eℋ1,1,3E>a/bBMS1β<μaIC13xo<a/b<Eℋ1,1,4E>a/bBMS1β<μaIC14xo=Eℋ1,1,5E>a/bBMS1β<μaIC15xo=a/bℋ1,2,1E>a/bBMS2β>μaIC21E<xoℋ1,2,2E>a/bBMS2β>μaC22a/b<xo<Eℋ1,2,3E>a/bBMS2β>μaIC23xo<a/b<Eℋ1,2,4E>a/bBMS2β>μaIC24xo=Eℋ2,2,5E>a/bBMS2β>μaIC25xo=a/bj=1,2 sign discriminates ordering between β and μa, and the composing index ki differentiates the position of the initial condition among the phase regions constituting a model type M1.Table 2
Basic phase arrangementsℋ2,j,ki determined by the PLF-PM and phase arrangements for model type M2 associated to the E<a/b ordering.
ℋi,j,kiM2BMSiBMSiICikiICikiℋ2,1,1E<a/bBMS1β<μaIC11E<xoℋ2,1,2E<a/bBMS1β<μaIC12a/b<xo<Eℋ2,1,3E<a/bBMS1β<μaIC13xo<a/b<Eℋ2,1,4E<a/bBMS1β<μaIC14xo=Eℋ2,1,5E<a/bBMS1β<μaIC15xo=a/bℋ2,2,1E<a/bBMS2β>μaIC21E<xoℋ2,2,2E<a/bBMS2β>μaC22a/b<xo<Eℋ2,2,3E<a/bBMS2β>μaIC23xo<a/b<Eℋ2,2,4E<a/bBMS2β>μaIC24xo=Eℋ2,2,5E<a/bBMS2β>μaIC25xo=a/bj=1,2 sign discriminates ordering between β and μa, and the composing index ki differentiates the position of the initial condition among the phase regions constituting a model type M2Then, for instance, the statementℋ2,2,1∧ICλi,j∧ICλ1,1≡1 signifies that we have a model type M1 with β>μa and that the initial condition satisfies a/b<xo, the λ1 equilibrium in region ℜ1x, that is, λ1>a/b, and finally that a/b<xo<λ1. In the Appendix, we present an extended qualitative study of the performance of the global trajectory xt under conjunctions ℋi,j,ki∧Oλi,j∧ICλi,j.Figure1(a) shows the performance of the PLF-PM composite trajectory xt for the ℋ1,1,1 configuration. Particularly, shown placement a/b<λ1<E derives from the auxiliary ordering E>μa2/βb. The global trajectory xt initiates in a region ℜ1x according to the x1t rule. Then, since for β<μa population size decreases, continuity of linked trajectory projects a time t1 such that xt reaches the E threshold, that is, xt1=E, and then xt enters into the region ℜ2x switching to the x2t rule. The settings imply the existence of a second time t2 so that xt touches the a/b frontier, thereby getting into the region ℜ3x and following according to the x3t rule. But, since β<μa, eventually xt vanishes. ℋ1,1,2 and ℋ1,1,3 configurations similarly tied to the β<μa ordering also drive population size to vanish Figure 1(b) pertains to the ℋ2,1,1 arrangement. Composite trajectory xt begins in the region ℜ1x according to the x1t rule, but β<μa, and then population size decreases, so at a time t1, it reaches the a/b threshold, that is, xt1=a/b, entering into the region ℜ4x and subsequently abiding by the rule x4t. The decreasing trend implies population size touching the E boundary at a second time t2, thereby placing within the region ℜ3x and then switching again, this time to the x4t rule. Afterward, population size keeps decreasing following an asymptotic trend to extinction. The faith of the PLF-PM composite trajectory xt for ℋ2,1,2 and ℋ2,1,3 is extinction equally. Figure 1(c) pertains to the ℋ3,1,1 conjunction corresponding to the xo>E ordering. Then, in the beginning, population size xt places in the region ℜ1x and decreases according to the x1t rule. Subsequently, it asymptotically approaches the value λ1<E, but on its trend, it hits the E boundary at time t1 switching to rule x3t holding in region ℜ3x. The condition β<μa drives population to extinction. For ℋ3,1,2 conjunction, the initial condition placement x0<E keeps the population size within region ℜ3x where it progresses to extinction according to the x3t trajectory.Figure 1
Examples of the construction of the PLF-PM composite trajectoryxt for ℋi,1,1 for i=1,2,3; these configurations include a BMS1:β<μa condition. (a) Composite trajectory xt for ℋ1,1,1. (b) Composite trajectory xt for ℋ2,1,1. (c) Composite trajectory xt for ℋ3,1,1.
(a)(b)(c)Figure 2
Examples of the construction of the PLF-PM composite trajectoryxt for ℋi,2,k for i=1,2, and k=1,2,3; referred configuration entitles the BMS2:β>μa condition. (a) The statement λ2>E∧ℋ1,2,1≡1 holds, λ1 is a stable equilibrium in the region ℜ1x. (b) Under the λ2<E∧ℋ1,2,2≡1 statement, λ2 places in region ℜ2x and becomes a stable equilibrium solution. (c) Whenever E<λ4<a/b∧ℋ2,2,2≡1, the λ4 equilibrium places in ℜ4x; since E<xo<a/b, xt increases, reaching the a/b boundary at time t1, and switches to rule x1t but since λ1>a/b, xt increases asymptotically towards λ1 staying within region ℜ1x.(d) In turn, if E<λ4<a/b∧ℋ2,2,2≡1, then λ4 stands for a stable equilibrium solution in the region ℜ4x.
(a)(b)(c)(d)Figure 3
Examples of the shape of the global trajectoryxt for ℋi,3,kiconfiguration. (a) The xt path for the ℋ1,3,1 conjunction population size xt starts in the region ℜ1x and approaches asymptotically the equilibrium solution λ2=a/b. (b) For ℋ2,3,2, population size xt begins in the region ℜ1x according to x1t and decreases and then switches to rule x4t to approach asymptotically λ4=E. (c) ℋ3,3,1 ordering sets xt initially within the region ℜ1x and then xt decreases to λ1=E. (d) For ℋ3,3,2, we have x0<E, and population size starts within the region ℜ3x and follows growth law x3t but since β=μa, it remains stationary at x0.
(a)(b)(c)(d)AssumeOλ2,1∧ℋ1,2,1≡1; then, the condition BMS1:β<μa holds; additionally, the λ2 equilibrium localizes above E, that is, Oλ2,1:λ2>E, and these orderings appear in conjunction with M1:E>a/b and IC11:E<xo. Particularly, λ2>E implies λ1>E so the initial condition placement E<xo<λ1sets the xt trajectory to increase according to the x1t rule, approaching the asymptotic limit λ1 and staying within region ℜ1x. On the other hand, whenever we have E<λ1<xo, population size xt decreases asymptotically towards λ1. This establishes λ1 as a stable equilibrium in region ℜ1x (see Figure 2(a)). Alternatively, under the statement Oλ2,2∧ℋ1,2,2≡1, since we have β>μa, the ordering λ2>a/b; then, since Oλ2,2 implies λ2<E, we could have a/b<xo<λ2<E or a/b<λ2<xo<E. If a/b<xo<λ2<E gets true, xt will behave according to the x2t law, and then it will decrease from xo and approach the asymptotic limit λ2.For a/b<λ2<xo<E, the population size xt still conforming to the x2t growth law will increase from xo and approach λ2.Thus, λ2 stands for a stable equilibrium in the region ℜ2x (see Figure 2(b)). Under the ℋ2,2,2 arrangement, population size xt begins in the region ℜ4x following the x4t rule with λ4>E. If E<x0<λ4<a/b,xt increases as it asymptotically approaches λ4. On the other hand, if E<λ4<x0<a/b, xt decreases while asymptotically approaching λ4. Therefore, the condition E<λ4<a/b∧ℋ2,2,2≡1 bears a stable equilibrium at λ4 (see Figure 2(c)). Whenever λ4>a/b, the λ4 equilibrium lies within region ℜ1x. Since E<xo<a/b, population size increases according to rule x4t approaching λ4 reaching the a/b boundary at time t1 switching to rule x1t, since λ1>a/b population size increases asymptotically towards λ1 staying within region ℜ1x (see Figure 2(d)). ℋ1,2,3 and ℋ2,2,3 bear xt increasing from E<xo and could approach asymptotically λ4 whenever this equilibrium solution lies within region ℜ4x or λ1 if λ4>a/b.Figure 4
Fit of the PLF-PPM to data of Armstrong [7] relates to the growth of an experimental population of planarian Dugesia tigrina. (a) Fitted trajectory lines. (b) Proportion 0.0260×Rx of resource abatement function.
(a)(b)Forℋ1,3,1, population size xt starts in the region ℜ1x and controls by x1t with equilibrium solution λ1 satisfying a/b<λ1<E<x0. Hence, xt decreases and attains the E threshold at a time t1, transferring to rule x2t and then asymptotically approaching the equilibrium solution λ2=a/b (see Figure 3(a)). For the arrangement ℋ1,3,2, the initial condition sets xt to acquire the x2t law in ℜ2x and then it approaches the equilibrium solution λ2. Correspondingly, the ℋ1,3,3 composite xt obeys the x3t rule in region ℜ3x, but since β=μa, it remains stationary at xo. These configurations suggest a maximum limiting effect of resource scarcity. For ℋ2,3,1, population size xt begins within the region ℜ1x and since λ1<a/b places in region ℜ4x, then x1t decreases, intersecting a/b at a time t1 where it switches to rule x4t, ultimately approaching the equilibrium solution λ4=E. A ℋ2,3,2 arrangement sets E<x0<a/b, and population size starts within the region ℜ4x and sets by x4t approaching the λ4=E equilibrium (see Figure 3(b)). For ℋ2,3,3, the conditions x0<E and β=μa set population size beginning within the region ℜ3x, following x3t that remains stationary at x0. For ℋ3,3,1, initial condition sets population size xt within the region ℜ1x, following x1t. Since β=μa and E=a/b implies λ1=E, then xt decreases to λ1 (see Figure 3(c)). For ℋ3,3,2, we have x0<E, and population size starts within the region ℜ3x and follows to growth law x3t but since β=μa, it remains stationary at x0 (see Figure 3(d)). Configurations in Figure 3 suggest a maximum limiting effect of resource scarcity on population growth; even a BMS1:β<μa condition does not perform.Figure 5
Fit of the PLF-PPM to data taken from Huisman [8] relates to the growth of an experimental population of the unicellular green algae, Chlorella vulgaris. (a) Fitted trajectory lines. (b) Proportion 0.0949×Rx of resource abatement function.
(a)(b)In summary, the PLF-PPM’s continuous-time global population size trajectory can stand an outstanding array of different performances. It could model population decline whenever the orderingβ<μa holds (see Figure 1). Alternatively, β>μa could induce conditional stability depending on the achieved parametric arrangement (e.g., Figure 2). Also, the statement β>μa could entail a purely growing regime for population size xt (e.g.E<λ4<a/b∧ℋ2,2,2≡1 shown in Figure 2(d). Furthermore, β=μa circumstance signposts a maximum effect of resource abatement on population growth (e.g., ℋ3,3,1 shown in Figure 3(c)). Moreover, for critical values of E threshold, this configuration could exacerbate Allee effects, and this is due to random influences which could promote vanishing of population size. Furthermore, the ratio a/b in the settings of the PLF-PPM is interpreted as an Allee threshold since xt<a/b implies axt>bx2t, so population size below a/b controls the mortality process as entailed by the Mx setup (cf. equation (21)). It is also pertinent to emphasize that for model type M1, even though population size places below the starvation threshold E, its dynamics could take place under an Allee effect regime. Correspondingly, for M2 type models, population could evolve in a regime combining both starvation and Allee effects. Consequently, by adapting the number of factors defining Nx and Mx in equation (17) so as to suit specific modeling aims, we could endure the PLF-PPM suitable predictive strength. In what follows, we explain the performance of the PLF-PPM model of equation (22) as an exploratory tool given different datasets.


## 3.1. Analytic Exploration

To study the varied forms of the global trajectoryxt derived from the PLF-PM, we begin by setting the array of possible phase configurations associated with the basic parameters a,b,E,β, and μ (cf. equations (22) and (23)). We depart from the ordering relationship for the external energy E and the ratio of mortality from population density proportions, a/b, that yields a first model type classification Mi,i=1,2,3. We have to consider models M1 associated to the inequality E>a/b, where the range of xt composes the regions: ℜ1x=a/b<E<x, ℜ2x=a/b<x<E, and ℜ3x=x<a/b<E (cf. equation (23)). For model type M2 where the relationship E>a/b holds, the range of xt will be divided into the regions ℜ1x,ℜ3x, and ℜ4x=x<a/b<E (cf. equation (23)). For model type M3 linking to the ordering E=a/b, the range of xt will be divided into the regions ℜ1x and ℜ3x. We can be aware that regions ℜ2x and ℜ4x are mutually excluding. Besides, it also requires classifying possible phase arrangements according to a Birth to Mortality Scaled Ordering BMSi,i=1,2,3, namely, BMS1: when β<μa, BMS2: if β>μa, and BMS3: whenever β=μa. Additionally, each model type Mi with i=1,2 determines possible placements of the initial condition, and each one associated to composing phase regions ℜix with i=1,2,3,4 determines an initial condition ordering ICikiwithi=1,2andki=1,2,…,5, namely, IC11: for (E<xo), IC12: whenever a/b<xo<E, IC13: when xo<a/b, IC14: if xo=E, and IC15: whenever xo=a/b. Correspondingly, IC21: for xo<E, IC22: whenever E<xo<a/b, IC23: when xo>a/b,IC24: if xo=E, and IC25: whenever xo>a/b. The possible phase arrangements typify according to a three-dimensional conjunction operator ℋi,j,ki=Mi∧BMSj∧ICiki for iand ki as specified. Table 1 summarizes phase arrangements for model type M1, and Table 2 summarizes those corresponding to model type M2. In addition to the ℋi,j,ki conjunctions, we must consider the positioning of the equilibrium solutions inside their associated regions. For this reason, in the presentation of results, we could use the Oλi,j indicator for i=1,2,4 designating equilibrium λi and j=1,2,4 describing region ℜix. Then, for instance, Oλ1,1 identifies the placement of equilibrium λ1 inside region ℜ1x; correspondingly, Oλ2,3 labels positioning of λ2 inside region ℜ3x and so on. Yet given Oλi,j, this relates to distinguishing how the initial condition is positioned relative to the equilibrium solution λi. Resulting orderings generically symbolize through ICλi,j for i=1,2,4 associated to the equilibrium solution λi and j=1,2 labeling the ordering of x0 relative to λi. This way, ICλi,1 stands for x0<λi and ICλi,2 stands for x0>λi. Besides described symbols, we could refer to direct inequalities that a given equilibrium satisfies relative to thresholds defining regions.Table 1
Basic phase arrangementsℋ1,j,ki determined by the PLF-PM and phase arrangements for model type M1 associated to the E>a/b ordering.
ℋi,j,kiM1BMSiBMSiICikiICikiℋ1,1,1E>a/bBMS1β<μaIC11E<xoℋ1,1,2E>a/bBMS1β<μaIC12a/b<xo<Eℋ1,1,3E>a/bBMS1β<μaIC13xo<a/b<Eℋ1,1,4E>a/bBMS1β<μaIC14xo=Eℋ1,1,5E>a/bBMS1β<μaIC15xo=a/bℋ1,2,1E>a/bBMS2β>μaIC21E<xoℋ1,2,2E>a/bBMS2β>μaC22a/b<xo<Eℋ1,2,3E>a/bBMS2β>μaIC23xo<a/b<Eℋ1,2,4E>a/bBMS2β>μaIC24xo=Eℋ2,2,5E>a/bBMS2β>μaIC25xo=a/bj=1,2 sign discriminates ordering between β and μa, and the composing index ki differentiates the position of the initial condition among the phase regions constituting a model type M1.Table 2
Basic phase arrangementsℋ2,j,ki determined by the PLF-PM and phase arrangements for model type M2 associated to the E<a/b ordering.
ℋi,j,kiM2BMSiBMSiICikiICikiℋ2,1,1E<a/bBMS1β<μaIC11E<xoℋ2,1,2E<a/bBMS1β<μaIC12a/b<xo<Eℋ2,1,3E<a/bBMS1β<μaIC13xo<a/b<Eℋ2,1,4E<a/bBMS1β<μaIC14xo=Eℋ2,1,5E<a/bBMS1β<μaIC15xo=a/bℋ2,2,1E<a/bBMS2β>μaIC21E<xoℋ2,2,2E<a/bBMS2β>μaC22a/b<xo<Eℋ2,2,3E<a/bBMS2β>μaIC23xo<a/b<Eℋ2,2,4E<a/bBMS2β>μaIC24xo=Eℋ2,2,5E<a/bBMS2β>μaIC25xo=a/bj=1,2 sign discriminates ordering between β and μa, and the composing index ki differentiates the position of the initial condition among the phase regions constituting a model type M2Then, for instance, the statementℋ2,2,1∧ICλi,j∧ICλ1,1≡1 signifies that we have a model type M1 with β>μa and that the initial condition satisfies a/b<xo, the λ1 equilibrium in region ℜ1x, that is, λ1>a/b, and finally that a/b<xo<λ1. In the Appendix, we present an extended qualitative study of the performance of the global trajectory xt under conjunctions ℋi,j,ki∧Oλi,j∧ICλi,j.Figure1(a) shows the performance of the PLF-PM composite trajectory xt for the ℋ1,1,1 configuration. Particularly, shown placement a/b<λ1<E derives from the auxiliary ordering E>μa2/βb. The global trajectory xt initiates in a region ℜ1x according to the x1t rule. Then, since for β<μa population size decreases, continuity of linked trajectory projects a time t1 such that xt reaches the E threshold, that is, xt1=E, and then xt enters into the region ℜ2x switching to the x2t rule. The settings imply the existence of a second time t2 so that xt touches the a/b frontier, thereby getting into the region ℜ3x and following according to the x3t rule. But, since β<μa, eventually xt vanishes. ℋ1,1,2 and ℋ1,1,3 configurations similarly tied to the β<μa ordering also drive population size to vanish Figure 1(b) pertains to the ℋ2,1,1 arrangement. Composite trajectory xt begins in the region ℜ1x according to the x1t rule, but β<μa, and then population size decreases, so at a time t1, it reaches the a/b threshold, that is, xt1=a/b, entering into the region ℜ4x and subsequently abiding by the rule x4t. The decreasing trend implies population size touching the E boundary at a second time t2, thereby placing within the region ℜ3x and then switching again, this time to the x4t rule. Afterward, population size keeps decreasing following an asymptotic trend to extinction. The faith of the PLF-PM composite trajectory xt for ℋ2,1,2 and ℋ2,1,3 is extinction equally. Figure 1(c) pertains to the ℋ3,1,1 conjunction corresponding to the xo>E ordering. Then, in the beginning, population size xt places in the region ℜ1x and decreases according to the x1t rule. Subsequently, it asymptotically approaches the value λ1<E, but on its trend, it hits the E boundary at time t1 switching to rule x3t holding in region ℜ3x. The condition β<μa drives population to extinction. For ℋ3,1,2 conjunction, the initial condition placement x0<E keeps the population size within region ℜ3x where it progresses to extinction according to the x3t trajectory.Figure 1
Examples of the construction of the PLF-PM composite trajectoryxt for ℋi,1,1 for i=1,2,3; these configurations include a BMS1:β<μa condition. (a) Composite trajectory xt for ℋ1,1,1. (b) Composite trajectory xt for ℋ2,1,1. (c) Composite trajectory xt for ℋ3,1,1.
(a)(b)(c)Figure 2
Examples of the construction of the PLF-PM composite trajectoryxt for ℋi,2,k for i=1,2, and k=1,2,3; referred configuration entitles the BMS2:β>μa condition. (a) The statement λ2>E∧ℋ1,2,1≡1 holds, λ1 is a stable equilibrium in the region ℜ1x. (b) Under the λ2<E∧ℋ1,2,2≡1 statement, λ2 places in region ℜ2x and becomes a stable equilibrium solution. (c) Whenever E<λ4<a/b∧ℋ2,2,2≡1, the λ4 equilibrium places in ℜ4x; since E<xo<a/b, xt increases, reaching the a/b boundary at time t1, and switches to rule x1t but since λ1>a/b, xt increases asymptotically towards λ1 staying within region ℜ1x.(d) In turn, if E<λ4<a/b∧ℋ2,2,2≡1, then λ4 stands for a stable equilibrium solution in the region ℜ4x.
(a)(b)(c)(d)Figure 3
Examples of the shape of the global trajectoryxt for ℋi,3,kiconfiguration. (a) The xt path for the ℋ1,3,1 conjunction population size xt starts in the region ℜ1x and approaches asymptotically the equilibrium solution λ2=a/b. (b) For ℋ2,3,2, population size xt begins in the region ℜ1x according to x1t and decreases and then switches to rule x4t to approach asymptotically λ4=E. (c) ℋ3,3,1 ordering sets xt initially within the region ℜ1x and then xt decreases to λ1=E. (d) For ℋ3,3,2, we have x0<E, and population size starts within the region ℜ3x and follows growth law x3t but since β=μa, it remains stationary at x0.
(a)(b)(c)(d)AssumeOλ2,1∧ℋ1,2,1≡1; then, the condition BMS1:β<μa holds; additionally, the λ2 equilibrium localizes above E, that is, Oλ2,1:λ2>E, and these orderings appear in conjunction with M1:E>a/b and IC11:E<xo. Particularly, λ2>E implies λ1>E so the initial condition placement E<xo<λ1sets the xt trajectory to increase according to the x1t rule, approaching the asymptotic limit λ1 and staying within region ℜ1x. On the other hand, whenever we have E<λ1<xo, population size xt decreases asymptotically towards λ1. This establishes λ1 as a stable equilibrium in region ℜ1x (see Figure 2(a)). Alternatively, under the statement Oλ2,2∧ℋ1,2,2≡1, since we have β>μa, the ordering λ2>a/b; then, since Oλ2,2 implies λ2<E, we could have a/b<xo<λ2<E or a/b<λ2<xo<E. If a/b<xo<λ2<E gets true, xt will behave according to the x2t law, and then it will decrease from xo and approach the asymptotic limit λ2.For a/b<λ2<xo<E, the population size xt still conforming to the x2t growth law will increase from xo and approach λ2.Thus, λ2 stands for a stable equilibrium in the region ℜ2x (see Figure 2(b)). Under the ℋ2,2,2 arrangement, population size xt begins in the region ℜ4x following the x4t rule with λ4>E. If E<x0<λ4<a/b,xt increases as it asymptotically approaches λ4. On the other hand, if E<λ4<x0<a/b, xt decreases while asymptotically approaching λ4. Therefore, the condition E<λ4<a/b∧ℋ2,2,2≡1 bears a stable equilibrium at λ4 (see Figure 2(c)). Whenever λ4>a/b, the λ4 equilibrium lies within region ℜ1x. Since E<xo<a/b, population size increases according to rule x4t approaching λ4 reaching the a/b boundary at time t1 switching to rule x1t, since λ1>a/b population size increases asymptotically towards λ1 staying within region ℜ1x (see Figure 2(d)). ℋ1,2,3 and ℋ2,2,3 bear xt increasing from E<xo and could approach asymptotically λ4 whenever this equilibrium solution lies within region ℜ4x or λ1 if λ4>a/b.Figure 4
Fit of the PLF-PPM to data of Armstrong [7] relates to the growth of an experimental population of planarian Dugesia tigrina. (a) Fitted trajectory lines. (b) Proportion 0.0260×Rx of resource abatement function.
(a)(b)Forℋ1,3,1, population size xt starts in the region ℜ1x and controls by x1t with equilibrium solution λ1 satisfying a/b<λ1<E<x0. Hence, xt decreases and attains the E threshold at a time t1, transferring to rule x2t and then asymptotically approaching the equilibrium solution λ2=a/b (see Figure 3(a)). For the arrangement ℋ1,3,2, the initial condition sets xt to acquire the x2t law in ℜ2x and then it approaches the equilibrium solution λ2. Correspondingly, the ℋ1,3,3 composite xt obeys the x3t rule in region ℜ3x, but since β=μa, it remains stationary at xo. These configurations suggest a maximum limiting effect of resource scarcity. For ℋ2,3,1, population size xt begins within the region ℜ1x and since λ1<a/b places in region ℜ4x, then x1t decreases, intersecting a/b at a time t1 where it switches to rule x4t, ultimately approaching the equilibrium solution λ4=E. A ℋ2,3,2 arrangement sets E<x0<a/b, and population size starts within the region ℜ4x and sets by x4t approaching the λ4=E equilibrium (see Figure 3(b)). For ℋ2,3,3, the conditions x0<E and β=μa set population size beginning within the region ℜ3x, following x3t that remains stationary at x0. For ℋ3,3,1, initial condition sets population size xt within the region ℜ1x, following x1t. Since β=μa and E=a/b implies λ1=E, then xt decreases to λ1 (see Figure 3(c)). For ℋ3,3,2, we have x0<E, and population size starts within the region ℜ3x and follows to growth law x3t but since β=μa, it remains stationary at x0 (see Figure 3(d)). Configurations in Figure 3 suggest a maximum limiting effect of resource scarcity on population growth; even a BMS1:β<μa condition does not perform.Figure 5
Fit of the PLF-PPM to data taken from Huisman [8] relates to the growth of an experimental population of the unicellular green algae, Chlorella vulgaris. (a) Fitted trajectory lines. (b) Proportion 0.0949×Rx of resource abatement function.
(a)(b)In summary, the PLF-PPM’s continuous-time global population size trajectory can stand an outstanding array of different performances. It could model population decline whenever the orderingβ<μa holds (see Figure 1). Alternatively, β>μa could induce conditional stability depending on the achieved parametric arrangement (e.g., Figure 2). Also, the statement β>μa could entail a purely growing regime for population size xt (e.g.E<λ4<a/b∧ℋ2,2,2≡1 shown in Figure 2(d). Furthermore, β=μa circumstance signposts a maximum effect of resource abatement on population growth (e.g., ℋ3,3,1 shown in Figure 3(c)). Moreover, for critical values of E threshold, this configuration could exacerbate Allee effects, and this is due to random influences which could promote vanishing of population size. Furthermore, the ratio a/b in the settings of the PLF-PPM is interpreted as an Allee threshold since xt<a/b implies axt>bx2t, so population size below a/b controls the mortality process as entailed by the Mx setup (cf. equation (21)). It is also pertinent to emphasize that for model type M1, even though population size places below the starvation threshold E, its dynamics could take place under an Allee effect regime. Correspondingly, for M2 type models, population could evolve in a regime combining both starvation and Allee effects. Consequently, by adapting the number of factors defining Nx and Mx in equation (17) so as to suit specific modeling aims, we could endure the PLF-PPM suitable predictive strength. In what follows, we explain the performance of the PLF-PPM model of equation (22) as an exploratory tool given different datasets.


## 4. Study Cases

In what follows, we explain the performance of the PLF-PPM model of equation (22) as an exploratory tool given different datasets. We address data of Armstrong [7] on the growth of populations of asexual Dugesia tigrina, data reported by Huisman [8] on developing experimental populations of the unicellular green algae Chlorella vulgaris, a study by Davidson [9] on a sheep population introduced in Tasmania, data by Pearl [10] on the growth of Drosophila melanogaster, and data reported by Hughes and Tanner [11] on the slow decline of an Agaricia agaricites population on Jamaican reefs. Fitted parameters, associated standard deviations, and concordance correlation coefficient values [12] are given in Table 3. All required PLF-PPM fits were achieved by using the Berkeley Madonna Software Version 8.3.18. Besides, acquiring the resource abatement function Rx relied on using equation (18) setting E=F0/2k for k=1, which assures getting a proxy for the maximum depletion rate.Table 3
Estimated values of initial population sizex9 and basic parameter a,b,β,μ,E, produced by fitting the PLF-PPM of equation (22) to the listed datasets.
DatasetFitted parametersxoAbBμEa/bλ1λ2λ4CCCArmstrong [7]34.313.590.090.030.00503077.4339.8847.410079.916.5e + 030.99200.060.380.060.000.0003625.18—————Huisman [8]0.75186.5339.710.680.00101744.692.12e + 316.512.1e + 30.99460.1149.6515.310.020.000284.56—————Pearl [10]3.961.890.080.250.0092628.3022.1838.27318.289.0e−30.99131.150.480.000.010.000137.48—————Davidson [9]85.731460.560.181.94e + 41974260.711.02e + 31.69e + 31.2e + 40.948458.1041.791.750.020.0000178.63—————Hughes and Tanner [11]213.421186.561.570.091.64e + 4376.55755.77180.74360.98752.050.93401.5153309.591.950.059.49e − 5104.61—————Records right below parameter estimates correspond to associated standard errors. Calculateda/b ratio, equilibria λ1,λ2,λ4, and concordance correlation coefficient (CCC) values are also displayed.Armstrong [7] maintained populations of asexual Dugesia tigrina, in an arrangement of finger bowls. Bowls contained 120 ccs of water each and were kept at a temperature of 25°C. Every other day, each population received 0.1 ccs of freshly killed brine shrimp. The bowls were cleaned at the end of the feeding period. The initial population size in every bowl amounted to 35 worms. Reproduction occurred only by transverse fission, with each worm dividing to produce a tail that developed into a new individual. For data assembly, a tail was any recent fission product not adequately developed to consume food. In the experiment, tails were added to the population’s artificially increasing reproductive efficiency. Because this method amounts to exogenous addition of biomass, we focused on data of population control named 1–0 without tails added to fit the PLF-PM. The plot of the CCC=0.9920 global trajectory is shown in Figure 4(a). Associated estimated parameter values are presented in Table 3. Acquired orderings were β>μa and xo<a/b<λ1<λ2<E<λ4 corresponding to a ℋ1,2,3 arrangement composing regions ℜ1x,ℜ2x, and ℜ3x. Population size places initially within the region ℜ3x and acquires the form x3t until it reaches the upper boundary at a/b at time t1=4.43 where it gets into the region ℜ2x switching to rule x2t, approaching the asymptotic limit λ2. Since dynamics mainly describe the region ℜ2x, mortality could perhaps be explained by cannibalism induced by high intraspecific competition. The author inferred about food scarcity for population levels near the equilibrium. Nevertheless, the PLF-PM identifies an energy threshold E well above the equilibrium level, which explains the shape of the resource depletion trajectory proxy shown in Figure 4(b). Therefore, population regulation is solely controlled by density-dependent mortality.Figure 6
Fit of the PLF-PPM to data of Davidson [9] relates to the growth of a sheep population introduced in Tasmania. (a) Fitted trajectory lines. (b) Resource abatement function R(x).
(a)(b)Huisman [8] studied the growth of experimental populations of the unicellular green algae Chlorella vulgaris. Populations were raised in continuous cultures incubated at 20°C and administered with plentiful nutrients and O2. Thus, growth limitation is exclusively associated with incident light energy. We handpicked data on population G. Figure 5(a) shows the global trajectory associated with the CCC=0.9946 fit of the PLF-PM. The values of estimated parameters are presented in Table 3. Acquired parametric orderings were β>μa and xo<a/b<λ2<E conforming to ℋ1,2,3 composing regions ℜ1x,ℜ2x, and ℜ3x. Population size places initially within the region ℜ3x and rules by x3t until it reaches the upper boundary at a/b at time t1=3.7 where it goes into the region ℜ2x switching to rule x2t, approaching the asymptotic limit λ2. The equilibrium value λ2 was practically attained at the end of the 18th day, with regulation by recruitment of new individuals and mortality dominated by individual properties. This photosynthetic organism synthesizes biomass based on a combination of factors like light, CO2, and nutrients. The threshold E could be conceived as the biomass size such that the aggregated life enduring factors become limiting factors. Because there were sufficient CO2 and nutrients in this case, the Ethreshold could be interpreted as the level of available light energy. Our analysis sustains the assumption in [8] that population growth in Chlorella vulgaris is limited exclusively by incident light. Therefore, the value of the Ethreshold estimated by the fit of the PLF-PPM determines a resource abatement function that suggests that the amount of incident light for the G-labeled population in Huisman’s experiment was not limiting (see Figure 5(b)); this is because the equivalent of E in biomass units lies well above the equilibrium level λ2.Figure 7
Fit of the PLF-PPM to data of Pearl [10] relates to the growth of an experimental population of the fruit fly Drosophila melanogaster. (a) Fitted trajectory lines. (b) Proportion 0.5066×Rx of resource abatement function Rx.
(a)(b)Davidson [9] studied a sheep population introduced in Tasmania. We acquired related proxy data from [13] representing the averages of the number of individuals taken on periods of 5 years. The recorded data span a total of 120 years. The first average corresponds to 1814–1819, and we arranged for 1819 to stand for year zero. The global trajectory acquired from the CCC=0.9484 fit of the PLF-PM is shown in Figure 6(a). The obtained parametric orderings were β>μa and xo<a/b<λ1<λ2<E<λ4, associated to a ℋ1,2,3 arrangement that composes regions ℜ1x,ℜ2x, and ℜ3x. Population size places initially within the region ℜ3x and progressed according to rule x3t. The fit points out that within the region ℜ3x, population growth was slow at the beginning, which suggests that reproduction was limited by a reduced number of individuals. Population size crossed the a/b threshold at the beginning of year 6 (1825). Population growth switched to rule x2t, associated with recruitment and mortality controlled by population size. Thus, fitted trajectory increased approaching the equilibrium λ2 placing inside ℜ2x (see Figure 6(a)). Besides, it can be learned from Figure 6(b) that in spite of surpassing two records and another two almost reaching the E threshold, their influence was not strong enough so as to change the damped oscillating trend about λ2 imposed by rule x2t. Indeed, the dominance of the dynamic arrangement endured by the x2t rule stayed over some time long enough to induce the trend to equilibrium suggested by the PLF-PPM fit.Figure 8
Fit of the PLF-PM to data of Hughes and Tanner [11] relates to the decline of Agaricia agaricites coral population on Jamaican reefs. (a) Fitted trajectory lines. (b) Resource abatement function Rx.
(a)(b)We now consider the construction of the global trajectoryxt associated with a fit of the PLF-PM to data of [10], taken from [14]. Pearl [10] maintained populations of Drosophila melanogaster in bottles using yeast as food and fitted a logistic curve to the associated data expressed as the counted number of adult flies in a particular day. Sang [15] criticized the procedure, quarreling that the experiments did not keep the yeast constant but as a growing population on its own. The plots associated with the PLF-PPM fit displayed in Figure 7(a) produced the parameter orderings β>μa and x0<a/b<E<λ2 that conform to a ℋ1,2,3 arrangement, that is, a model type M1 with BMS2: β>μa and IC13: xo<a/b. Phase arrangement includes regions ℜ1x,ℜ2x, and ℜ3x. At the beginning of the growth process, xt stays in the region ℜ3x but the statement BMS2: β>μa≡1, forces associating a rule x3t to increase, and then it reaches the a/b boundary at t1=10, thereby entering the region ℜ2x and then progressing asymptotically to λ2 following rule x2t. The fit shows that the fly population is entirely controlled by density-dependent mortality. In contrast to Sang [15], the PLF-PPM fit demonstrates that energy did not play a decisive role E>λ2; also, see Figure 7(b)), so we can close that variations in feeding energy did not determine the dynamics since otherwise x1t would compose the global trajectory xt. The application of the FDM sustains the assumption that in Pearl’s experiment, energy does not provide a criterion to define a carrying capacity. Furthermore, the equilibrium λ2 is not expressed in terms of E but solely as a function of the scaling parameters β,μ, and b.Figure 9
The PLF-PPM composite trajectoryxt for the ℋ1,1,1 configuration. (a) a/b<βE/μb<E derived from the auxiliary ordering E>μa2/βb. This panel shows that the global trajectory xt reaches the E threshold at a time t1. (b) Population size keeps decreasing, so it touches the a/b frontier. (c) Since β<μa, the population becomes extinct.
(a)(b)(c)Hughes and Tanner [11] reported a slow decline of the coral population on Jamaican reefs over 16 years. We fitted the PLF-PM to data on Agaricia agaricites colonies recorded during 1977–1993. The authors established that local extinction is explained by increased mortality rates plus the adverse effects of two hurricanes, Allen (in 1980) and Gilberto (in 1988). Also relevant in explaining the decline was the impossibility of recovery because of a recruitment failure. Figure 8(a) shows the global trajectory associated with the fit CCC=0.9340 of the PLF-PM. The values of the parameter estimates are given in Table 3. Acquired parametric orderings were β<μa and xo<λ2<E<λ1<a/b conforming to a ℋ2,1,3 arrangement and composing regions ℜ1x,ℜ2x, and ℜ3x. Population size places initially within the region ℜ3x and rules by x3t. Food was not limiting. Perhaps perturbation set its level under the upper boundary of the ℜ3x region where Allee effects dominate mortality. Then, possibly a recruitment failure impeded compensation to the abatement effects linked to the β<μa ordering. Then, our analysis corroborates the assumption of Hughes and Tanner [11]. In any event, recruitment failure could not associate with resource reduction (see Figure 8(b)).Figure 10
The PLF composite trajectoryxt for the ℋ1,1,2 configuration. (a) The global trajectory xt decreases from xo until it reaches the a/b threshold at time t1. Since β<μa, population size keeps decreasing until it becomes extinct. (b) Extinction as determined by the ℋ1,1,3 arrangement.
(a)(b)The PLF-PPM displayed high reproducibility strength in all performed fits, as it derives from the high CCC values included in Table3. Moreover, parallel plots of resource abatement explain that fitted values of initial resource availability thresholds E explain the occurrence of recorded equilibrium levels. Results point towards the consistency of the PLF-PPM as an exploratory tool. Nevertheless, the direct fitting procedure endured by the Berkeley Madonna software relied upon apportioning initial values of parameter estimates. Necessarily, this brought about high sensibility associated to local minimum problems at the non-linear acquisition of final estimates. Detected inconveniences suggest revision aimed at adapting techniques that could lessen the experienced parameter estimation burden.


### 4.1. Outline of the Discrete-Time Setup of the PLF-PPM

Although this paper focuses on addressing the analysis of the PLF-PPM arrangement trajectories on continuous time, we consider it pertinent at this setting to outline the associated discrete-time form. For this reason, we letfixn denote the discrete-time form of the phase model fix as given by equation (22) for i=1,2,3,4. Also, let ℜixn be the corresponding range of application of the rule fixn. Then, if χℜ1xn stands for the characteristic function of ℜ1xn, we get(25)xn+1=xn+∑i=14χℜixnfixn,for n=0,1,2,…, with x0 being placed on one of the regions ℜixn. Then,(26)f1xn=βE−μbxn2ℜ1xn=xn|xn>E∧xn>abf2xn=βxn−μbxn2ℜ2xn=xn|xn<E∧xn>abf3xn=β−μaxnℜ3xn=xn|xn<E∧xn<abf4xn=βE−μaxnℜ4xn=xn|xn>E∧xn<ab.The parametrizationr=1+β and s=μb and then the change of variable r/szn=xn will bring the logistic map as a submodel composing the seek for discrete-time formulation. Concomitantly, equation (25) transforms into(27)zn+1=∑1=1mχℜiznθizn,where(28)θ1zn=s1−rEr+zn1−rzn,ℜ1zn=zn|zn>sEr∧zn>sarbθ2zn=rzn1−znℜ2zn=zn|zn<sEr∧zn>sarbθ3zn=r−sabznℜ3zn=zn|zn<sEr∧zn<sarbθ4zn=sr−1Er+1−sabznℜ4zn=zn|zn>sEr∧zn<sarb.Equivalently,(29)zn+1=s1−rEr+zn1−rzn,ℜ1zn=zn>sEr∧zn>sarb,rzn1−zn,ℜ2zn=zn<sEr∧zn>sarb,r−sabzn,ℜ3zn=zn<sEr∧zn<sarb,sr−1Er+1−sabzn,ℜ4zn=zn>sEr∧zn<sarb.For the sake of conciseness, we avert a detailed study of the performance of the discrete-time form of the PLF-PPM to an upcoming second part of this work. In the meantime, the parametrization of equation (25) given by equation (29) bears insight since it suggests that given the condition on the parameter r that renders upset of chaos in the logistic map submodel θ2zn, once the resulting orbit crosses the upper boundary of the region ℜ2zn, dynamics will be set by the θ1zn submodel so the system will escape the chaotic regime and approach an equilibrium point in the region ℜ1zn or the θ1zn orbit could get back again to the chaotic regime determined by the logistic map θ2zn. Alternatively, the θ2zn orbit could reach the lower boundary of the region ℜ2zn then switching to the θ3zn dynamics ruling in ℜ3zn. Then depending in the paramter odering zn could decrease steadily and eventually vanish, or else wander aound in region R3zn before getting back to the chaotic regime.


## 4.1. Outline of the Discrete-Time Setup of the PLF-PPM

Although this paper focuses on addressing the analysis of the PLF-PPM arrangement trajectories on continuous time, we consider it pertinent at this setting to outline the associated discrete-time form. For this reason, we letfixn denote the discrete-time form of the phase model fix as given by equation (22) for i=1,2,3,4. Also, let ℜixn be the corresponding range of application of the rule fixn. Then, if χℜ1xn stands for the characteristic function of ℜ1xn, we get(25)xn+1=xn+∑i=14χℜixnfixn,for n=0,1,2,…, with x0 being placed on one of the regions ℜixn. Then,(26)f1xn=βE−μbxn2ℜ1xn=xn|xn>E∧xn>abf2xn=βxn−μbxn2ℜ2xn=xn|xn<E∧xn>abf3xn=β−μaxnℜ3xn=xn|xn<E∧xn<abf4xn=βE−μaxnℜ4xn=xn|xn>E∧xn<ab.The parametrizationr=1+β and s=μb and then the change of variable r/szn=xn will bring the logistic map as a submodel composing the seek for discrete-time formulation. Concomitantly, equation (25) transforms into(27)zn+1=∑1=1mχℜiznθizn,where(28)θ1zn=s1−rEr+zn1−rzn,ℜ1zn=zn|zn>sEr∧zn>sarbθ2zn=rzn1−znℜ2zn=zn|zn<sEr∧zn>sarbθ3zn=r−sabznℜ3zn=zn|zn<sEr∧zn<sarbθ4zn=sr−1Er+1−sabznℜ4zn=zn|zn>sEr∧zn<sarb.Equivalently,(29)zn+1=s1−rEr+zn1−rzn,ℜ1zn=zn>sEr∧zn>sarb,rzn1−zn,ℜ2zn=zn<sEr∧zn>sarb,r−sabzn,ℜ3zn=zn<sEr∧zn<sarb,sr−1Er+1−sabzn,ℜ4zn=zn>sEr∧zn<sarb.For the sake of conciseness, we avert a detailed study of the performance of the discrete-time form of the PLF-PPM to an upcoming second part of this work. In the meantime, the parametrization of equation (25) given by equation (29) bears insight since it suggests that given the condition on the parameter r that renders upset of chaos in the logistic map submodel θ2zn, once the resulting orbit crosses the upper boundary of the region ℜ2zn, dynamics will be set by the θ1zn submodel so the system will escape the chaotic regime and approach an equilibrium point in the region ℜ1zn or the θ1zn orbit could get back again to the chaotic regime determined by the logistic map θ2zn. Alternatively, the θ2zn orbit could reach the lower boundary of the region ℜ2zn then switching to the θ3zn dynamics ruling in ℜ3zn. Then depending in the paramter odering zn could decrease steadily and eventually vanish, or else wander aound in region R3zn before getting back to the chaotic regime.


## 5. Discussion

Regulation by extreme value is observed in a plentiful of biological processes covering varied scales. It manifests on cellular structures such as mitochondria, where the maxima or minima of a periodical chemical reaction rule the formation of observed patterns [16]. Also, in the ecological settings, variables such as physical stress due to high or low temperatures, salinity, soil water content, wind velocities, and long or short exposures to air express better through extreme values than standard measures of central tendency [17, 18]. The notion of biological control by extreme values dates back to Justus Von Liebig, who established the law of the minimum. It states that an organism’s growth rate is regulated by the nutrient present in the minimum [5]. Generally, any factor that slows down potential growth in an ecosystem is described as a limiting factor. Acknowledgment of lower-upper tolerance limits for a process drove generalization of the law of the minimum into the law of the tolerance of Shelford [6].Furthermore, the interaction of limiting factors can indirectly influence the effect of other factors not in themselves limiting. This paradigm, known as the principle of limiting factors, has proven to be very useful in studying whole or parts of ecosystems [13, 19]. The joint Liebig–Shelford paradigm adopted here bears that population growth control occurs by balancing birth and mortality processes, the first determined by factor inducing the minimum offspring and the second by forcing the maximum number of deaths among individuals. This notion sustains the formal piecewise setup expressed by equation (22) and refers to PLF-PPM. Such a paradigm could allow a piecewise account of population dynamics composing growth phases delimited by density dependence cooperation and competition [20–23], thresholds of starvation, and critical density or extinction [20, 24]. The present PLF-PPM conceives growth phases controlled by limiting factors acting over domains bounded by population size thresholds. The first version of a PLF-PPM addressed a predator-prey model where Liebig’s law governs the natality process of the prey population [1]. Echavarría and Gomez [25] and Echavarría et al. [26] extended these ideas to formulate models in which the mortality rate is set by the maximum value of factors that depend on population size. Montiel-Arzate et al. [3] adopted the referred Liebig–Shelford principle of limiting factors to formulate a functionally diverse population growth model. The approach in [3] relied on a parametrization of the birth process that rendered qualitative exploration burdensome, explaining why it is missing. Besides, Montiel-Arzate et al. [3] did not include a specific parametrization of mortality due to Allee effects which we incorporate here represented by the term axt as described by equation (21). Present settings allowed exhaustive qualitative exploration of the PLF-PPM global trajectory xt presented in the Appendix. Even considering only population size and an external resource as limiting factors, the present formulation ensures a set of 20 basic orderings, each one associated to a different distribution of phases, which conform to an outstanding array of varying operation modes for the PLF-PPM. Such a flexible structure endures a sound interpretative strength displayed in the fitting results of the addressed study cases. Particularly, present settings allow visualization of abatement of resources controlling the birth process.Nevertheless, a direct fitting procedure that relies on estimates’ initial values brought high sensibility associated with local minimum problems at the non-linear acquisition of final values. These inconveniences suggest revision aimed to adapt a maximization of reproducibility strength criterion and simulation techniques that could lessen the parameter estimation burden experienced. Another point concerns adaptation of a symbolic manipulation code aimed at automatically exploring phase arrangement and their implications for stability. Yet another essential issue left untouched here concerns exploring the outlined discrete-time version of the PLF-PPM of equation (25) and its alternate form (29). In the meantime, the parametrization of equation (29) already suggests that even conditions on the parameter r rendering upset of chaos in the logistic map submodel θ2zn once the resulting orbit crosses the upper boundary of the regionℜ2zn dynamics will be set by the θ1zn submodel, so the system will escape chaos before the orbit gets back again to the chaotic regime determined by the logistic map θ2zn. Given the present results, no doubt that the suggested adaptation of codes to enhance both qualitative exploration and parameter fitting procedures will strengthen the applicability of offered PLF-PPM. Meanwhile, its present form already provides an outstanding research tool through which critical ecological parameters can be identified and be meaningfully interpreted.


## 6. Conclusions

This paper elucidates a logical deductive approach to establishing a piecewise structured model to interpret the growth of a single species population. The offered PLF-PPM derives as a logical consequence of what we coined as a principle of limiting factors for population growth, a paradigm adapted by merging Liebig’s law of the minimum and the tolerance law of Shelford. The formal approach explains by the extreme value characterization of the birthNxt and mortality Mxt processes in equation (15). Conceived forms allow the selection of the factors that play the decisive role in controlling population growth at a given time. The flexibility from a piecewise structure imbues the PLF-PPM with a noticeable reproducibility strength. Considered study cases briefly elaborate on this. The PLF-PPM bears as well an outstanding interpretative advantage. Such a feature infuses our construct with the capability of identifying different growth phases. These associate with regions of the dominance of crucial factors determining population dynamics such as a starvation regime, population size viability, or regimes of high population density effects. We only contemplated one external factor, E, and population size itself to describe Nxt and Mxt. Concomitantly, the number of possible phase configurations was over ten times the number of parameters involved, which is the reason why ensuing parametrization proved to be seemingly adequate given addressed datasets. Nevertheless, as the number of explaining factors increases, more parameters are required. Then, the tied analytical exploration becomes complicated. Then particularly, acquiring suitable complexity could lead to local minimum difficulties connected to non-linear estimation. In this order of ideas, a quest for efficient parameter estimation methods that enhance the practical advantages derived from the PLF-PPM seems necessary. Such an endeavor concerns the research aims to pursue in a further paper. Another relevant research subject is pending, which is exhaustively exploring the discrete-time PLF-PPM construct outlined by equation (25).



---

*Source: 5623783-2021-12-27.xml*