Technique to derive discrete population models with delayed growth

Abstract

We provide a procedure for deriving discrete population models for the size of the adult population at the beginning of each breeding cycle and assume only adult individuals reproduce. This derivation technique includes delay to account for the number of breeding cycles that a newborn individual remains immature and does not contribute to reproduction. These models include a survival probability (during the delay period) for the immature individuals, since these individuals have to survive to reach maturity and become members of, what we consider, the adult population. We discuss properties of this class of discrete delay population models and show that there is a critical delay threshold. The population goes extinct if the delay exceeds this threshold. We apply this derivation procedure to obtain two models, a Beverton–Holt adult model and a Ricker adult model and discuss the global dynamics of both models.

Keywords:

• Delay difference equations
• global stability
• persistence
• critical delay threshold
• Beverton–Holt
• Ricker

1. Introduction

Difference equation models of populations are a well-established modelling tool to describe populations and are used as an underlying modelling framework to study the dynamics of species and derive sustainable control strategies such as harvest regulations [8, 13, 14, 33, 34]. Formulating discrete population models requires a clear understanding of the length of the underlying time interval between consecutive time points. Once the interpretation of time units is determined, the processes that affect the population have to be described with respect to these time units. For example, if a species reproduces once a year in a particular month and the time unit is assumed to be one month, then the reproductive event only occurs every 12th time point and should not be present in the equation for every time point. In most applications, the time points refer to consecutive years and a reproductive event is assumed to occur once in that year. Thus, the length between time points matches the time between consecutive breeding cycles. In that case, in classical formulations of the form $X_{t+1}=F\left(X_t\right)$, where F is a growth function, it is implicitly assumed that all individuals at time $t,X_t$, contribute to reproduction. Thus, newborn individuals at time $t$ are assumed to contribute to the reproduction and therefore the growth of the population after only one time step, i.e. newborn individuals are assumed to reach (sexual) maturity within one breeding cycle.

However, several, especially long-lived species reach (sexual) maturity after several breeding cycles (and therefore several years) after hatching. For example, Barramundi fish, also known as Asian Bass, are only able to reproduce after age 4–6 [10]. Similarly, Leatherback turtles reach maturity on average at age 24 [7] and Bowhead whales reach maturity after about 23 years [25]. Clearly, these are only a few examples of a wide variety of species that require several breeding cycles to reach maturity. In these cases, delay equations have been argued to be more appropriate [9, 12, 24, 27, 28].

In this manuscript, we propose a derivation technique to model the subgroup of mature individuals. Based on the general form $X_{t+1}=F\left(X_t\right)$, we first express the growth function $F\left(X_t\right)$ as a combination of surviving (sexually) mature individuals and individuals that join the group of mature individuals for the first time at time $t+1$. Thus, we separate the growth function with respect to cohorts that contribute to the growth of the mature population. This separation allows us to implement the model assumption that individuals require several breeding cycles until reaching maturity. Since the immature individuals have to survive until the age they reach maturity, we incorporate a survival probability for each cohort of immature individuals (see Section 2).

The resulting discrete population models differ from existing formulations of discrete delay population models. See for example [1, 4, 9, 11, 19, 21, 27, 30]. Although this is only a small collection of existing delay recurrences in discrete population modelling, it highlights the different derivation strategies. For example, the model in [19] used a commonly exploited method to formulate discrete population models, namely the application of the Euler discretization scheme to a continuous (delay) population model. Although commonly applied, this derivation method can result in models with negative solutions, making it biologically questionable. Another popular class of delay difference equations in population modelling assumes a delay in the fitness function, defined by the ratio of solutions evaluated at consecutive time points. This takes the form $X_{t+1}=X_{t}F\left(X_{t-\tau }\right)$, also known as the Levin-May model [18] (see for example [21]). In fact, the well-known Pielou model [27] is an example of such a formulation. Alternatively, one can derive discrete delay population models in a bottom-up style based on an age-structured population model, as was done in [4]. We will see that this bottom-up derivation corresponds to a delay recurrence that can account for a delay in maturation.

A related and widely accepted class of delay population models was proposed in [9] to account for delay in maturity, observed in Antarctic fin whales. The proposed delay recurrence takes the general form $X_{t+1}=\alpha X_t+F\left(X_{t-\tau }\right)$ and was studied by several authors. For example, in [11], the dynamics of such models was studied and several examples of the recruitment function F were provided. This proposed delay structure considers delay in the reproductive function, while assuming no delay in the survival. In this work, we build on this idea of introducing delay in the reproductive function. This was also the strategy behind the derivation in [1, 30]. However, in [1], only conditions for the local asymptotic stability of the equilibrium are provided and the effect of delay is not studied. Also, in [30], only a specific growth model, the Beverton–Holt survival function, was considered.

In this work, we instead formulate a class of delay recurrences for general survival functions, therefore extending the idea in [30]. As well, even though in both formulations, it is assumed that the delay is only in the reproduction, the assumptions that underpin the derivation of the class of models in this paper differ from the assumptions in [30], resulting in a different model structure. While the assumptions in [30] result in a recurrence of the form $X_{t+1}=F_1\left(X_t\right)+F_2(X_t,X_{t-\tau })$, our assumptions in this manuscript lead to $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)$. The class of delay difference equations that we derive in Section 2 has realistic model properties. For example, solutions remain non-negative for non-negative initial population values. Our analysis of the general model, in Section 3, shows that there is a model specific delay threshold. If the delay exceeds this delay threshold, then the population goes extinct. This is a biologically realistic property, as the length of time to reach maturity increases, the probability that hatched individuals survive until maturity naturally declines. This in turn, reduces the overall contribution to actual growth.

In Section 4, we derive two specific population models using the technique introduced. In the first model, we derive a Beverton–Holt adult population model, where both the mature individuals and the cohort of immature individuals are exposed to their cohort-specific Beverton–Holt survival function. We prove the global dynamics by showing that if the delay remains below the critical delay threshold, then the population converges to a unique positive equilibrium. For the second-model, the Ricker-adult population model, we assume that mature individuals have a Ricker survival function and each cohort of immature individuals is exposed to a Beverton–Holt survival function. We show that if the delay is below the model specific delay threshold, the population persists and discuss the local asymptotic stability of the unique positive equilibrium. We extend this result and provide sufficient conditions for the global asymptotic stability of the positive equilibrium, but also show that it is possible for this model to have an attracting period two solution when the positive equilibrium is unstable.

2. Derivation of a discrete delay growth model

Many discrete models of population growth in the literature, in which $X_t$ can be interpreted as the density of individuals in the population at the beginning of the time interval between $t$ and $t+1$, are of the form (1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) for $t\in {\mathbb{N}}_0=\{0,1,2,3,\dots \}$. Expanding the right-hand side of (1(1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) ) yields $X_{t+1}=p\left(X_t\right)X_t+p\left(X_t\right)g\left(X_t\right)X_t=F_1\left(X_t\right)+F_2\left(X_t\right).$One can interpret $X_t$ as density of adult individuals that are sexually mature at the beginning of the time interval $t$ to $t+1$ (i.e. the end of the cycle from $t-1$ to $t$). Here, $p\left(X_t\right)\in (0,1)$ is the survival probability of those individuals that are alive at time $t$ (i.e. the beginning of the interval from $t$ to $t+1$) and survive to time $t+1$ (i.e. to the end of the interval from $t$ to $t+1$). Let $F_1\left(X_t\right):=p\left(X_t\right)X_t$.  Then $F_1\left(X_t\right)$ represents the number of individuals in the population alive and sexually mature at the end of the interval from $t$−1 to $t$ that survive the time interval from $t$ to $t+1$. The term $F_2\left(X_t\right):=p\left(X_t\right)g\left(X_t\right)X_t$ describes the density of individuals that were born at the beginning of the breeding cycle from  t to $t+1$ that reach maturity at the end of this cycle (i.e. at time $t+1$). In this formulation, eggs laid at the beginning of the interval from $t$ to $t+1$ are assumed to produce sexually mature individuals that are ready to reproduce within one time interval (i.e. at the beginning of the time interval from $t+1$ to $t+2$).

One could easily incorporate the assumption that the survival probability of the individuals just born at the beginning of the interval from $t$ to $t+1$ is different from the survival probability of already mature individuals. Then, model (1(1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) ) can be replaced by the slightly more general description (2) $X_{t+1}=p\left(X_t\right)X_t+\hat{p}\left(X_t\right)g\left(X_t\right)X_t,$(2) where $\hat{p}\in (0,1)$ represents the survival probability during one time interval of newborn individuals.

The well-established discrete single species model (3) $X_{t+1}=\frac{\rho K{X}_t}{K+(\rho -1)X_t},$(3) also known as Beverton–Holt model was introduced in [5] for K>0, $\rho >1$. This model is of form (1(1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) ) with (4) $p\left(X_t\right)=\frac{1}{1+(\rho -1)\frac{X_t}{K}}=\hat{p}\left(X_t\right),g\left(X_t\right)=\rho -1>0.$(4) Another popular discrete single species model, the Ricker model, (5) $X_{t+1}=X_t{e}^{r-d-c{X}_t},$(5) for $r$, $d$, $c>0$, introduced in [29], is also of the form in (1(1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) ) with (6) $p\left(X_t\right)=e^{-d-c{X}_t}=\hat{p}\left(X_t\right),g\left(X_t\right)=(e^r-1)>0.$(6) An implicit assumption in (1(1) $X_{t+1}=p\left(X_t\right)(1+g(X_t\left)\right)X_t,$(1) ) and (2(2) $X_{t+1}=p\left(X_t\right)X_t+\hat{p}\left(X_t\right)g\left(X_t\right)X_t,$(2) ) is that newborn individuals join the population and are able to procreate at the beginning of the breeding cycle immediately following the one in which they were born. However, this assumption is inaccurate for many species as maturity is often reached only after several breeding cycles. Thus, the reproduction is subject to a delay that takes into consideration that it can require more than one breeding cycle for newborns to reach maturity and be able to reproduce. Accounting for such a delay will also require a modification of the survival probability of the newborns because individuals that are born at the beginning of the breeding season $t-\tau$, where $\tau \in {\mathbb{N}}_0$, can only reproduce for the first time at the beginning of the breeding season $t+1$ if they survive this delay period. This implies that it takes individuals $\tau +1$ breeding cycles to reach maturity and contribute to reproduction.

In order to refine the general model (2(2) $X_{t+1}=p\left(X_t\right)X_t+\hat{p}\left(X_t\right)g\left(X_t\right)X_t,$(2) ) to take this delay into account, we let $X_t$ denote the density of sexually mature individuals that can lay eggs at the beginning of the season from time $t$ to $t+1$. We let $\tau +1$ denote the number of breeding seasons from when the eggs are laid until the time they reach sexual maturity. Since the incubation time of the eggs is usually insignificant compared to the length of time it takes hatchlings to reach maturity, it will be ignored. Assuming also that the immature individuals only interact with other immature individuals of the same cohort, the model with delay that recognizes the need for hatchlings to survive to maturity before they can reproduce is given by (7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7)

where $h(\tau ,X_{t-\tau })=\tilde{p}(\tau ,X_{t-\tau })g\left(X_{t-\tau }\right)$ with $\tau +1$ initial conditions, $X_i\ge 0$, $i\in \{0,1,\dots ,\tau \}$.

Again, the first term, $F_1\left(X_t\right)$ contains the already mature individuals at time $t$ that survive to time $t+1$ and $F_2\left(X_{t-\tau }\right)=h(\tau ,X_{t-\tau })X_{t-\tau }$ captures the cohort of previously immature individuals that reach, for the first time, sexual maturity at time $t+1$. As in the model without delay, $p\left(X_t\right)\in (0,1)$ is the survival probability of adults over one breeding season and so $p\left(X_t\right)X_t$ denotes the density of adults that were alive at the end of the cycle from $t$−1 to $t$ that survive to the end of the cycle from $t$ to $t+1$. The term, $\tilde{p}(\tau ,X_{t-\tau })\in (0,1)$ is the survival probability of immature individuals born (and hatched) in the beginning of the time interval $t$ that survive to time $t+1$, when they first join the adult population and hence are part of the breeding population. Therefore, $\tilde{p}(\tau ,X_{t-\tau })g\left(X_{t-\tau }\right)X_{t-\tau }$ represents the density of eggs laid at time $t-\tau$ that hatch and survive to become adults at time $t+1$ and can reproduce for the first time at time $t+1$ (i.e. the beginning of the cycle from $t+1$ to $t+2$).

Model (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) could, of course, be generalized by considering that the fraction of immature individuals that survive the delay period depends not only on the same cohort but rather on all of the individuals (both mature and immature individuals). In this case, $\tilde{p}(\tau ,X_{t-\tau })$ would be replaced by $\tilde{p}(\tau ,X_{t-\tau },X_{t-\tau +1},\dots ,X_t)$. Nevertheless, the specific assumption in (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) can be understood as an approximation where the competition with individuals of the same cohort is stronger than with individuals of other cohorts. This implicit assumption has also been exploited by works in [3, 9, 24, 27].

It is worth noting that (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) can be derived from a stage-structured population model of the form ${\overrightarrow{Y}}_{t+1}=B_t{\overrightarrow{Y}}_t$, for ${\overrightarrow{Y}}_t=(A_t,J_t^{\left(1\right)},\dots ,J_t^{\left(\tau \right)})$, where $({\overrightarrow{Y}}_t{)}_1=A_t$ represents the mature individuals and $({\overrightarrow{Y}}_t{)}_i=J_t^{(i-1)}$ for $i=2,\dots ,\tau +1$, where $J^{\left(i\right)}$ represents the ith age-class of immature individuals with $J^{\left(1\right)}$ representing the newborns that survive the breeding cycle they are born in. Assuming that the survival of immature individuals in each age class to the next age-class is cohort dependent, that is, the survival from class $J^{\left(i\right)}$ to $J^{(i+1)}$ is dependent on $J^{\left(i\right)}$, the density dependent Leslie Matrix is then given by the $(\tau +1)\times (\tau +1)$ matrix (8) $B_t=\left[\begin{array}{cccccc}p\left(A_t\right)& 0& 0& \dots & 0& q_{\tau }\left(J_t^{\left(\tau \right)}\right)\\ q_0\left(A_t\right)& 0& 0& \dots & 0& 0\\ 0& q_1\left(J_t^{\left(1\right)}\right)& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & q_{\tau -1}\left(J_t^{(\tau -1)}\right)& 0\end{array}\right],$(8) where $p\left(z\right)\in (0,1)$ denotes the fraction of already mature individuals z that survive from time $t$ to $t+1$ and for $i>0$, $q_i\left(w\right)\in (0,1)$ denotes the fraction of immature individuals w in state i, that survive to the next breeding cycle. The function $q_0\left(A\right)$ represents the newborn individuals that survive the time interval that they are born in.

Note that there is density dependence in all the non-zero terms in $B_t$, and $B_t$ differs from density dependent Leslie matrices that consider only the dependence on the total population, as in [2, 4, 15].

Focusing only on the mature adult population, this age-structured population model becomes, after recursive application, (9) $\begin{array}{ll}{A}_{t+1}& =p\left(A_t\right)A_t+q_{\tau }\left(J_t^{\left(\tau \right)}\right)J_t^{\left(\tau \right)}\\ & =p\left(A_t\right)A_t+q_{\tau }\left(q_{\tau -1}\right(J_{t-1}^{(\tau -1)}\left)J_{t-1}^{(\tau -1)}\right)q_{\tau -1}\left(J_{t-1}^{(\tau -1)}\right)J_{t-1}^{(\tau -1)}\\ & =⋮\\ & =p\left(A_t\right)A_t+h(\tau ,A_{t-\tau })A_{t-\tau },\end{array}$(9) for an appropriate function $h(\tau ,A_{t-\tau })$. However, by modelling the mature individuals and therefore considering (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ), we are able to exploit classical methods to investigate the impact of the length of the delay on the survival of the population (i.e. the number of immature age-classes). Instead, in an age-structured population model, this analysis would require the comparison of discrete systems with different system length, i.e. a sequence of matrices with increasing dimension.

Throughout, we consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) with $g\in C^1(\mathbb{R},{\mathbb{R}}^{+})$, $p\in C^1(\mathbb{R},(0,1\left)\right)$, and $\tilde{p}\in C^1({\mathbb{N}}_0\times \mathbb{R},(0,1\left)\right)$. In Section 3, we investigate this model with the additional assumptions that are all satisfied for the models discussed in Section 4:

(H1)	$p\left(u\right)$ is non-increasing in $uR$, i.e. $p_u:=\frac{\mathrm{d}p\left(u\right)}{\mathrm{d}u}\le 0$.
(H2)	$\tilde{p}(\tau ,u)$ is (strictly) decreasing (to zero) in $\tau$ such that ${\tilde{p}}_{\tau }:=\frac{\mathrm{\partial }\tilde{p}(\tau ,u)}{\mathrm{\partial }\tau }<0$ and $\underset{\tau \to \mathrm{\infty }}{lim}\tilde{p}(\tau ,u)=0$, and non-increasing in $uR$, i.e. ${\tilde{p}}_u:=\frac{\mathrm{\partial }\tilde{p}(\tau ,u)}{\mathrm{\partial }u}\le 0$.
(H3)	For fixed $\tau$, there exists a finite $u_{crit}\left(\tau \right)$ such that $\frac{\mathrm{\partial }h(\tau ,u)}{\mathrm{\partial }u}\ge 0$ for $u\le u_{crit}\left(\tau \right)$, and $\frac{\mathrm{\partial }h(\tau ,u)}{\mathrm{\partial }u}<0$ for $u_{crit}\left(\tau \right)\le u<\mathrm{\infty }$ and $\underset{\tau \to \mathrm{\infty }}{lim}h(\tau ,u_{crit}(\tau \left)\right)=0$.

The model assumptions (H1)–(H2) imply that there is no positive density effect on the survival probabilities. Thus, density independent survival probabilities or survival that is impacted by intra-specific competition are examples that satisfy $p_u,{\tilde{p}}_u\le 0$. Furthermore, it is reasonable to assume that with a longer delay period, the probability of surviving the delay period decreases, that is, ${\tilde{p}}_{\tau }<0$. Lastly, (H3) assumes that the per-capita growth function of the newly joining individuals is eventually decreasing for sufficiently large population densities. We will see that this condition is satisfied for popular examples, such as the Beverton–Holt and the Ricker model, where $h(\tau ,X)$ decreases for $X\ge 0$ (i.e. $X_{crit}<0$).

3. Model properties

In this section, we consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) with $\tau \in {\mathbb{N}}_0$ and non-negative initial conditions.

We introduce the following notation.

• Denote the set of ordered non-negative fixed points of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) by ${\mathcal{X}}^{\ast }=\{X_0^{\ast },X_1^{\ast },\dots \}$ with $X_i^{\ast }\le X_{i+1}^{\ast }$ for $i\in {\mathbb{N}}_0$. We assume that ${\mathcal{X}}^{\ast }$ is finite.

• We set $X_{max}^{\ast }=max\{X_0^{\ast },X_1^{\ast },\dots \}$ and assume $X_{max}^{\ast }<\mathrm{\infty }$.

• If (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) is assumed to satisfy (H3), we denote $X_{crit}$ to be the critical threshold replacing $u_{crit}$ in (H3).

• If $1-p\left(0\right)<h(0,0)$, then we let $\overline{\tau }>0$ be the unique value such that $1-p\left(0\right)=h(\overline{\tau },0)$.

Note that since $X=0$ is a fixed point, $X_0^{\ast }=0$.

The structure of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) immediately implies that solutions for non-negative initial conditions remain non-negative.

Proposition 3.1

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). If $X_i=0$ for all $i\in \{0,1,\dots ,\tau \}$, then $X_t=0$ for all $t\ge 0$. Otherwise, if $X_i>0$ for at least one $i\in \{0,1,\dots ,\tau \}$, then there exists $T>0$ such that $X_t>0$ for all $t\ge T$.

Proposition 3.2

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H2) hold and $h(\tau ,X)X$ is bounded, that is, $\underset{X\ge 0}{sup}h(\tau ,X)X<\mathrm{\infty }$. Solutions remain bounded for all $t\ge 0$.

Proof.

First, note that since $h(\tau ,X)X$ is bounded, there exists $M>0$ such that $h(\tau ,X)X<M$ for all $X>0$ and there exists $N>M$ such that $p\left(0\right)+M{N}^{-1}<1$.

Suppose the solution is not bounded. Then there exists a strictly increasing sub-sequence, $\left\{X_{t_i}\right\}$, such that $X_k<X_{t_i}$ for all $k<t_i$ and an $n$ such that $X_{t_i}>N$ for all $i>n$. Then, for all $i>n$, $\begin{array}{rl}{X}_{t_i}& =p\left(X_{t_i-1}\right)X_{t_i-1}+h(\tau ,X_{t_i-1-\tau })X_{t_i-1-\tau }<p\left(0\right)X_{t_i-1}+M{N}^{-1}N\\ & \le \left\{p\right(0)+M{N}^{-1}\}max\{X_{t_i-1},N\}<max\{X_{t_i-1},N\}<X_{t_i},\end{array}$resulting in a contradiction.

Proposition 3.3

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H2), $h(\tau ,X)X$ is bounded, and $X_0^{\ast }$ is either a repellor or a saddle with stable manifold that does not intersect the interior of the first quadrant. Assume $X_i>0$ for at least one $i\in \{0,1,\dots ,\tau \}$. There exists $ϵ>0$ and $T\ge 0$ such that $X_t\ge ϵ$ for all $t\ge T$.

Proof.

Let $X_t$ be any solution with non-negative and non-trivial (not all zero) initial conditions. By Proposition 3.1, there exists $T\ge 0$ such that $X_t>0$ for all $t\ge T$. Suppose, for the sake of contradiction that $\underset{t\to \mathrm{\infty }}{lim inf}{X}_t=0$. Then, there exists a sub-sequence, $\left\{X_{t_n}\right\}$ with $\underset{n\to \mathrm{\infty }}{lim}{X}_{t_n}=0$.

We show recursively $\tau$ times, passing to a sub-sequence each time if necessary, that this implies that for $n$ sufficiently large, $\underset{n\to \mathrm{\infty }}{lim sup}{X}_{t_n-k}=0$, $k=0,1,2,\dots ,\tau$ and hence $\underset{t\to \mathrm{\infty }}{lim}{X}_t=0$. Since, by assumption, the stable manifold does not intersect the interior of the first quadrant, this gives a contradiction.

In the first step, suppose for the original sub-sequence $\left\{X_{t_n}\right\}$ the sequence $\left\{X_{t_n-1}\right\}$ has $\underset{n\to \mathrm{\infty }}{lim sup}{X}_{t_n-1}>0$. Then, passing to a sub-sequence $\left\{X_{t_{n_1}}\right\}$ of $\left\{X_{t_n}\right\}$, if necessary, there exists a ${\delta }_1>0$ and $N_{{\delta }_1}>0$ such that $X_{t_{n_1}-1}>{\delta }_1$ for all $n_1\ge N_{{\delta }_1}$. Since $\underset{n\to \mathrm{\infty }}{lim}{X}_{t_n}=0$, we can choose ${ϵ}_1<p\left(\overline{X}\right){\delta }_1$ and $N_{{ϵ}_1}$ sufficiently large so that $X_{t_{n_1}}<{ϵ}_1$ for all $n_1>N_{{ϵ}_1}$. Here, $\overline{X}$ is the upper bound on solutions, given by Proposition 3.2 so that $X_t\le \overline{X}$ for all $t\ge 0$. Thus, for $N_1>max\{N_{{\delta }_1},N_{{ϵ}_1}\}$, $\begin{array}{rl}{X}_{t_{n_1}}& =p\left(X_{t_{n_1}-1}\right)X_{t_{n_1}-1}+h(\tau ,X_{t_{n_1}-1-\tau })X_{t_{n_1}-1-\tau }>p\left(X_{t_{n_1}-1}\right)X_{t_{n_1}-1}\ge p\left(\overline{X}\right)X_{t_{n_1}-1}\\ & \ge p\left(\overline{X}\right){\delta }_1.\end{array}$This contradicts $X_{t_{n_1}}<{ϵ}_1$ for $n_1>N_1$. Therefore, $\underset{n\to \mathrm{\infty }}{lim}{X}_{t_n}=0$ and $\underset{n\to \mathrm{\infty }}{lim}{X}_{t_n-1}=0$. Applying this argument recursively, $\tau$ times, and choosing $0<\delta <min\left\{{\delta }_i\right\},$ $0<ϵ<min\left\{{ϵ}_i\right\}$, and $i=1,2,\dots ,\tau$, it follows that $\underset{n\to \mathrm{\infty }}{lim}{X}_{t_n-k}=0,k=0,1,2\dots ,\tau$ and hence $\underset{t\to \mathrm{\infty }}{lim}{X}_t=0$, contradicting the assumption that the stable manifold of $X_0^{\ast }$ does not intersect the interior of the first quadrant.

Proposition 3.4

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H3) hold. Let $X_{crit}=X_{crit}\left(\tau \right)$.

(a)	If $h(\tau ,X_{crit})<1-p\left(X_{crit}\right)$, then $X_{max}^{\ast }<X_{crit}$. If $h(\tau ,X_{crit})=1-p\left(X_{crit}\right)$, then $X_{max}^{\ast }=X_{crit}$. Otherwise, if $h(\tau ,X_{crit})>1-p\left(X_{crit}\right)$, then $X_i^{\ast }<X_{crit}<X_{max}^{\ast }$ for all $X_i^{\ast }\in {\mathcal{X}}^{\ast }\mathrm{\setminus }\left\{X_{max}^{\ast }\right\}$.
(b)	If $X_{crit}\le 0$ and $1-p\left(0\right)\ge h(\tau ,0)$, then no positive equilibrium exists and ${\mathcal{X}}^{\ast }=\left\{X_0^{\ast }\right\}$.
(c)	If $X_{crit}\le 0$ and $h(\tau ,0)>1-p\left(0\right)$, then there exists a unique positive equilibrium $X_1^{\ast }$ and ${\mathcal{X}}^{\ast }=\{X_0^{\ast },X_1^{\ast }\}$.

Proof.

Clearly, $X_0^{\ast }$ is an equilibrium of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Let $X^{\ast }>0$ be an equilibrium of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Then, $1=p\left(X^{\ast }\right)+h(\tau ,X^{\ast })$, i.e. $1-p\left(X^{\ast }\right)=h(\tau ,X^{\ast })$. That is, $X^{\ast }$ is determined by the intersection of the functions $q_1\left(X\right)=1-p\left(X\right)$ and $q_2\left(X\right)=h(\tau ,X)$. Since $q_1\left(X\right)$ is increasing by (H2), and $q_2\left(X\right)$ is decreasing for $X>X_{crit}$ by (H3), $q_1\left(X\right)$ and $q_2\left(X\right)$ can only intersect for $X>X_{crit}$ if $q_1\left(X_{crit}\right)<q_2\left(X_{crit}\right)$. This completes the claim in a). Since $X_{crit}\le 0$ in b) and c), the proof is complete.

Theorem 3.5

Existence of a critical delay threshold for the stability of the extinction equilibrium

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H3). There exists a critical delay threshold ${\tau }_c$ such that if $\tau >{\tau }_c$, then ${\mathcal{X}}^{\ast }=\left\{X_0^{\ast }\right\}$ and $\underset{t\to \mathrm{\infty }}{lim}X\left(t\right)=0$.

Proof.

By (H3), there exists $X_{crit}\ge 0$ such that $\underset{X\ge 0}{max}h(\tau ,X)\le h(\tau ,X_{crit})$.

Together with (H1)–(H2), we have $\begin{array}{rl}{X}_{t+1}& =p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau }\le p\left(X_t\right)X_t+h(\tau ,X_{crit})X_{t-\tau }\\ & \le p\left(0\right)X_t+\tilde{p}(\tau ,0)g\left(X_{crit}\right)X_{t-\tau }\le \left(p\right(0)+\tilde{p}(\tau ,0\left)g\right(X_{crit}\left)\right)max\{X_t,X_{t-\tau }\}.\end{array}$By (H2), there exists $\hat{\tau }$ such that $p\left(0\right)+\tilde{p}(\hat{\tau },0)g\left(X_{crit}\right)<1$. Therefore, for $\tau \ge \hat{\tau }$, $X_{t+1}\le \left(p\right(0)+\tilde{p}(\tau ,0\left)g\right(X_{crit}\left)\right)max\{X_t,X_{t-\tau }\}=bmax\{X_t,X_{t-\tau }\}$for $b:=\left(p\right(0)+\tilde{p}(\tau ,0\left)g\right(X_{crit}\left)\right)<1$. Thus, taking ${\tau }_c=\hat{\tau }$, it follows, by [20, Theorem 2], that $X_0^{\ast }$ is globally asymptotically stable for solutions with non-negative initial conditions.

In Section 4, we provide examples for which (H1)–(H3) are satisfied and Theorem 3.5 applies.

Remark 3.6

If (H1)–(H3) hold, then the critical delay threshold ${\tau }_c$ is bounded above by $\overline{\tau }$, where $\overline{\tau }$ solves $p\left(0\right)+\tilde{p}(\overline{\tau },0)g\left(X_{crit}\right)=1$. Furthermore, if $p\left(0\right)+\tilde{p}(0,0)g\left(0\right)>1$, then $\overline{\tau }>0$.

Lemma 3.7

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H3) hold and that $0<X_{max}^{\ast }$ is finite. Then, $X_{max}^{\ast }=X_{max}^{\ast }\left(\tau \right)$ is decreasing in $\tau$.

Proof.

For notational convenience, let $x\left(\tau \right)$ represent $X_{max}^{\ast }\left(\tau \right)$. Then, (10) $1-p\left(x\right(\tau \left)\right)=h(\tau ,x(\tau \left)\right)$(10) and since we assume that (H1)–(H3) hold, $p^{\prime }\left(x\right(\tau \left)\right)+\frac{\mathrm{\partial }h(\tau ,x(\tau \left)\right)}{\mathrm{\partial }x}<0$.

Furthermore, taking the derivative with respect to $\tau$ on both sides of (10(10) $1-p\left(x\right(\tau \left)\right)=h(\tau ,x(\tau \left)\right)$(10) ) yields, $-p^{\prime }\left(x\right(\tau \left)\right)x^{\prime }\left(\tau \right)=\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }\tau }+\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }x}{x}^{\prime }\left(\tau \right).$After rearranging, we have $x^{\prime }\left(\tau \right)=\frac{-\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }\tau }}{\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }x}+p^{\prime }\left(x\right)}.$Note that $\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }\tau }=g\left(x\right)\frac{\mathrm{\partial }\tilde{p}(\tau ,x)}{\mathrm{\partial }\tau }<0,$and since $p^{\prime }\left(x\right)<0$, $\frac{\mathrm{\partial }h(\tau ,x)}{\mathrm{\partial }x}+p^{\prime }\left(x\right)<0$. Thus, $x^{\prime }\left(\tau \right)<0$, completing the proof.

Corollary 3.8

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H3). Let $g\left(X\right)=g>0$. If $1-p\left(0\right)<h(0,0)$, then the critical delay threshold ${\tau }_c=\overline{\tau }>0$, where $1-p\left(0\right)=h(\overline{\tau },0)$. If $\tau <{\tau }_c$, then there exists a unique positive equilibrium $X_{max}^{\ast }\left(\tau \right)$ that is decreasing in $\tau$. If $\tau >{\tau }_c$, the only non-negative equilibrium is $X_0^{\ast }$ and it is globally asymptotically stable.

In particular, if $1-p\left(0\right)\ge h(0,0)$, then $X_0^{\ast }$ is the only non-negative equilibrium for any $\tau >0$ and it is globally asymptotically stable.

Proof.

Since $g\left(X\right)=g$ and we assume (H3) holds, $X_{crit}<0$. If $1-p\left(0\right)<h(\tau ,0)$, by Proposition 3.4(c), there exists a unique positive equilibrium $X_{max}^{\ast }$. By Lemma 3.7, $X_{max}^{\ast }$ is decreasing in $\tau$. If instead $1-p\left(0\right)>h(\tau ,0)$, by Proposition 3.4(b), $X_{max}^{\ast }=X_0^{\ast }$. Since $h(\tau ,X)$ is decreasing in $\tau$ for fixed $X$, $X_{t+1}\le \left(p\right(X_t)+h(\tau ,X_{t-\tau }\left)\right)max\{X_t,X_{t-\tau }\}\le \left(p\right(0)+h(\tau ,0\left)\right)max\{X_t,X_{t-\tau }\}$and since $p\left(0\right)+h(\tau ,0)<1$, we have, by [20, Theorem 2], that $X_0^{\ast }$ is globally asymptotically stable.

To continue the discussion of the stability of the equilibria, we linearize (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) by rewriting (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) as a system of difference equations of length $\tau +1$ using $\overrightarrow{W}=(W^{\left(1\right)},W^{\left(2\right)},\dots ,W^{(\tau +1)})=(X_t,X_{t-1},\dots ,X_{t-\tau })$. Then, the Jacobian is of the form (11) $J=\left[\begin{array}{cc}a{\overrightarrow{e}}_1^T& b\\ I& 0\end{array}\right],$(11) where $I\in {\mathbb{R}}^{\tau \times \tau }$ is the identity matrix, ${\overrightarrow{e}}_1^T=(1,0,\dots ,0)\in {\mathbb{R}}^{\tau }$ and (12) $a=p\left(W^{\left(1\right)}\right)+p^{\prime }\left(W^{\left(1\right)}\right)W^{\left(1\right)},b=h(\tau ,W^{(\tau +1)})+W^{(\tau +1)}\frac{\mathrm{\partial }h(\tau ,v)}{\mathrm{\partial }v}{|}_{v=W^{(\tau +1)}}.$(12) The characteristic equation of (11(11) $J=\left[\begin{array}{cc}a{\overrightarrow{e}}_1^T& b\\ I& 0\end{array}\right],$(11) ) at an equilibrium $X$ is given by (13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) for $a,b$, defined in (12(12) $a=p\left(W^{\left(1\right)}\right)+p^{\prime }\left(W^{\left(1\right)}\right)W^{\left(1\right)},b=h(\tau ,W^{(\tau +1)})+W^{(\tau +1)}\frac{\mathrm{\partial }h(\tau ,v)}{\mathrm{\partial }v}{|}_{v=W^{(\tau +1)}}.$(12) ), evaluated at $W^{\left(i\right)}=X$ for $1\le i\le \tau +1$.

Proposition 3.9

Consider (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ). Assume (H1)–(H3) hold. If there exists ${\tau }_0$ such that $p\left(0\right)+h({\tau }_0,0)<1$, then $X_0^{\ast }$ is locally asymptotically stable for all $\tau \ge {\tau }_0$.

Proof.

By (13(13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) ), the characteristic equation at $X_0^{\ast }$ is given by ${\lambda }^{\tau +1}-p\left(0\right){\lambda }^{\tau }-h(\tau ,0)=0.$The assumption $p\left(0\right)+h({\tau }_0,0)<1$ is a sufficient condition for $X_0^{\ast }$ to be locally asymptotically stable as shown in [9]. Since $\frac{\mathrm{\partial }h(\tau ,0)}{\mathrm{\partial }\tau }=g\left(0\right)\frac{\mathrm{\partial }\tilde{p}(\tau ,0)}{\mathrm{\partial }\tau }<0$, $p\left(0\right)+h(\tau ,0)<1$ for all $\tau \ge {\tau }_0$. Hence, the result follows.

The previous result implies that there is a threshold for the stability of the trivial equilibrium with respect to $\tau$. That is, if $X_0^{\ast }\left({\tau }_0\right)$ is stable, then $X_0^{\ast }\left(\tau \right)$ remains stable for all $\tau \ge {\tau }_0$. It is not hard to show that if $X_0^{\ast }\left(0\right)$ is unstable, then there exists ${\tau }_0>0$ such that $X_0^{\ast }\left(\tau \right)$ is stable for all $\tau \ge {\tau }_0$ and unstable for $\tau <{\tau }_0$. The analogous statement for $X_{max}^{\ast }>0$ does not hold. That is, there does not have to be a unique delay threshold such that a positive equilibrium is stable for all  $\tau$ less than that threshold.

Consider for example (14) $X_{t+1}=\frac{X_t}{1+\sum _{i=0}^3{a}_i{X}_t^i}+\frac{r{X}_{t-\tau }}{1+\sum _{i=0}^3{c}_i\left(\tau \right)X_{t-\tau }^i}$(14) with $a_0=a_1=a_3=1,a_2=\frac{3}{2},c_0=\frac{1}{4}(1+\tau ),c_1=c_2=c_3=1,g\left(X_t\right)=r=55.$Then, (14(14) $X_{t+1}=\frac{X_t}{1+\sum _{i=0}^3{a}_i{X}_t^i}+\frac{r{X}_{t-\tau }}{1+\sum _{i=0}^3{c}_i\left(\tau \right)X_{t-\tau }^i}$(14) ) satisfies (H1)–(H3) with ${\tau }_c=435$. From the orbit diagram in Figure 1, the positive equilibrium is attractive for $50\le \tau \le {\tau }_c=435$. Figure 1 demonstrates that there is a positive equilibrium that is stable for relatively large values of $\tau$ i.e. $50<\tau <{\tau }_c$, but is unstable for relatively small values of $\tau$ i.e. $0<\tau <50$. Thus, in this example, increasing the delay in the model can stabilize an otherwise unstable positive equilibrium, although increasing it too much eliminates the positive equilibrium.

Figure 1. Orbit diagram for (14(14) $X_{t+1}=\frac{X_t}{1+\sum _{i=0}^3{a}_i{X}_t^i}+\frac{r{X}_{t-\tau }}{1+\sum _{i=0}^3{c}_i\left(\tau \right)X_{t-\tau }^i}$(14) ) created with Matlab [23]. (a) We plotted the last 1000 values of a total of 51000 iterations for randomly chosen positive initial conditions. $X_{\mathrm{\infty }}$ denotes the limiting behaviour of the simulated orbit. For small delay, there does not appear to be a stable equilibrium. Increasing the delay sufficiently, stabilizes a positive equilibrium. However, as predicted by Theorem 3.5, increasing the delay beyond ${\tau }_c$ results in extinction. (b) We plotted the last 100 values of a total of ${10}^7$ iterations for $\tau \in \{35,36,\dots ,55\}$ to demonstrate the type of dynamics that can occur when the positive equilibrium loses stability.

It is also worth noting that the local asymptotic stability of $X_0^{\ast }$ does not necessarily imply its global asymptotic stability. As an example, consider (15) $p\left(X_t\right)=\frac{1}{1.1+.5X_t^2},\tilde{p}(\tau ,X_{t-\tau })=\frac{1}{1.1+\frac{1}{2}{X}_{t-\tau }^2+2\tau },g\left(X_{t-\tau }\right)=10X_{t-\tau }.$(15) Then, (H1)–(H3) are satisfied and $h(1,X)=\tilde{p}(1,X)g\left(X\right)$ is increasing for $X<X_{crit}\approx 2.49$, and decreasing for $X>X_{crit}$. For $\tau =0$, (15(15) $p\left(X_t\right)=\frac{1}{1.1+.5X_t^2},\tilde{p}(\tau ,X_{t-\tau })=\frac{1}{1.1+\frac{1}{2}{X}_{t-\tau }^2+2\tau },g\left(X_{t-\tau }\right)=10X_{t-\tau }.$(15) ) gives the non-delayed population model $X_{t+1}=\frac{1+10X_t}{1.1+.5X_t^2}$ that was introduced in [31] as a population model with an Allee effect. For any $\tau \ge 0$, $X_0^{\ast }$ is locally asymptotically stable by [9], since $a=0.909091$ and $b=0$ in (13(13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) ). However, $X_0^{\ast }$ is not always globally asymptotically stable, since for example, for $\tau =1$, there exists two positive equilibria, $X_1^{\ast }\approx 0.0283$ (unstable) and $X_2^{\ast }\approx 19.789$ (locally asymptotically stable).

4. Applications

In this section, we apply the model derivation procedure introduced in Section 2 to two popular discrete growth models, the Beverton–Holt model and the Ricker model.

4.1. Beverton–Holt adult population model

Following the model derivation technique introduced in [31], we assume that the population at time $t+1$ is given by a fraction $f\left(t\right)=\frac{1+g\left(t\right)}{1+q\left(t\right)}$, where $g\left(t\right)$ describes the growth contribution and $q\left(t\right)$ the decline contributions. Thus, $X_{t+1}=\frac{1+g\left(t\right)}{1+q\left(t\right)}{X}_t.$We assume that the growth is determined by a constant rate $r>0$, and the decline rate is due to both the death rate $d>0$ and intraspecific competition with factor $c>0$. Then the single species model without delay is given by (16) $X_{t+1}=\frac{1+r}{1+d+c{X}_t}{X}_t.$(16) Choosing $\rho =\frac{1+r}{1+d}$ and $K=\frac{r-d}{c}>0$ transforms (16(16) $X_{t+1}=\frac{1+r}{1+d+c{X}_t}{X}_t.$(16) ) into the classical Beverton–Holt model $X_{t+1}=\frac{\rho K}{K+(\rho -1)X_t}{X}_t$ where $\rho >1$.

Expanding the terms in (16(16) $X_{t+1}=\frac{1+r}{1+d+c{X}_t}{X}_t.$(16) ) to obtain (17) $X_{t+1}=\frac{1+r}{1+d+c{X}_t}{X}_t=\frac{1}{1+d+c{X}_t}{X}_t+\frac{1}{1+d+c{X}_t}r{X}_t,$(17) we can identify the survival probability of mature individuals from time $t$ to $t+1$ as (18) $p\left(X_t\right)=\frac{1}{1+d+c{X}_t}$(18) and the growth contribution as $g\left(X_t\right)=r>0$.

Considering a delay in the growth contribution in the sense of Section 2, we immediately have (19) $X_{t+1}=\frac{1}{1+d+c{X}_t}{X}_t+\tilde{p}(\tau ,X_{t-\tau })r{X}_{t-\tau }.$(19) To determine the survival probability of a cohort of immature individuals from time $t-\tau$ to time $t+1$, we make the following assumptions:

1. The death rate of each cohort of immature individuals is $D>0$.

2. Each cohort of immature individuals compete only within their own cohort until the cohort reaches maturity. The competition factor is assumed to be $C>0$.

Based on the assumptions in (i) and (ii), the survival probability of each cohort of immature individuals within one time interval, say $t-\tau$ to $t-\tau +1$, is (20) $s_{t-\tau +1,t-\tau }=\frac{1}{1+D+C{w}_{t-\tau }}$(20) where $w_{t-\tau }$ is the density of immature individuals of the same cohort at time $t-\tau$. To determine the survival probability of this cohort of immature individuals from time $t-\tau$ to $t+1$, we solve the recurrence (21) $w_{t+1}=\frac{1}{1+D+C{w}_t}{w}_t.$(21) This is in the form of a Beverton–Holt equation and has the solution $w_t=\frac{D{w}_{t_0}}{D(1+D{)}^{t-t_0}+((1+D{)}^{t-t_0}-1)C{w}_{t_0}}.$Setting $t_0=t-\tau$, we obtain the fraction of immature individuals of the cohort born at the beginning of the time interval $(t-\tau ,t-\tau +1)$ that survive to time $t+1$: $w_{t+1}=\frac{D{w}_{t-\tau }}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)C{w}_{t-\tau }}.$Since, at time $t-\tau$, the cohort consisted of $g\left(X_{t-\tau }\right)X_{t-\tau }$ individuals, the initial condition is $w_{t-\tau }=g\left(X_{t-\tau }\right)X_{t-\tau }=r{X}_{t-\tau }$. Thus, the immature individuals of the cohort produced at the beginning of the time interval $(t-\tau ,t-\tau +1)$ that join the mature individuals for the first time at time $t+1$ are therefore $\tilde{p}(\tau ,X_{t-\tau })g\left(X_{t-\tau }\right)X_{t-\tau }=\frac{Dr{X}_{t-\tau }}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$Thus, (22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) Substituting the expression (22(22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) ) into (19(19) $X_{t+1}=\frac{1}{1+d+c{X}_t}{X}_t+\tilde{p}(\tau ,X_{t-\tau })r{X}_{t-\tau }.$(19) ) and setting $\text{β}:=(1+D{)}^{\tau +1}$ yields the delayed population model (23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) with positive model parameters $d,c,C$, and $r$, and initial conditions $X_i\ge 0$ for $i\in \{0,1,\dots ,\tau \}$. This delayed Beverton–Holt adult model (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) differs from existing delayed Beverton–Holt models such as the ones in [6, 27, 30].

Note that, based on the above model assumptions, the corresponding age-structured population model is of the form ${\overrightarrow{Y}}_{t+1}=B_t{\overrightarrow{Y}}_t$, for $B_t$ in (8(8) $B_t=\left[\begin{array}{cccccc}p\left(A_t\right)& 0& 0& \dots & 0& q_{\tau }\left(J_t^{\left(\tau \right)}\right)\\ q_0\left(A_t\right)& 0& 0& \dots & 0& 0\\ 0& q_1\left(J_t^{\left(1\right)}\right)& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & q_{\tau -1}\left(J_t^{(\tau -1)}\right)& 0\end{array}\right],$(8) ) with (24) $\begin{array}{rl}p\left(A\right)& \stackrel{\left(18\right)}{=}\frac{1}{1+d+cA},q_0\left(A\right)=\frac{r}{1+D+CrA},\\ q_i\left(w\right)& \stackrel{\left(20\right)}{=}\frac{1}{1+D+Cw},\text{for}1\le i\le \tau .\end{array}$(24) As in (8(8) $B_t=\left[\begin{array}{cccccc}p\left(A_t\right)& 0& 0& \dots & 0& q_{\tau }\left(J_t^{\left(\tau \right)}\right)\\ q_0\left(A_t\right)& 0& 0& \dots & 0& 0\\ 0& q_1\left(J_t^{\left(1\right)}\right)& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & q_{\tau -1}\left(J_t^{(\tau -1)}\right)& 0\end{array}\right],$(8) ), $p\left(A\right)$ represents the survival probability of mature individuals and $q_i\left(w\right)$, for $1\le i\le \tau$ represents the survival probability of immature individuals w from time $t$ to $t+1$. The function $q_0\left(A\right)$ represents the newborns that survive the time interval they are born in and corresponds, by the above derivation, to the first iteration of (21(21) $w_{t+1}=\frac{1}{1+D+C{w}_t}{w}_t.$(21) ) with initial condition $rA$, since this represents the density of newborn individuals.

Thus, $B_t$ in (8(8) $B_t=\left[\begin{array}{cccccc}p\left(A_t\right)& 0& 0& \dots & 0& q_{\tau }\left(J_t^{\left(\tau \right)}\right)\\ q_0\left(A_t\right)& 0& 0& \dots & 0& 0\\ 0& q_1\left(J_t^{\left(1\right)}\right)& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & q_{\tau -1}\left(J_t^{(\tau -1)}\right)& 0\end{array}\right],$(8) ) is $B_t=\left[\begin{array}{cccccc}\frac{1}{1+d+c{A}_t}& 0& 0& \dots & 0& \frac{1}{1+D+C{J}_t^{\left(\tau \right)}}\\ \frac{r}{1+D+rC{A}_t}& 0& 0& \dots & 0& 0\\ 0& \frac{1}{1+D+C{J}_t^{\left(1\right)}}& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & \frac{1}{1+D+C{J}_t^{(\tau -1)}}& 0\end{array}\right].$One can verify, that repeated substitution as suggested in (9(9) $\begin{array}{ll}{A}_{t+1}& =p\left(A_t\right)A_t+q_{\tau }\left(J_t^{\left(\tau \right)}\right)J_t^{\left(\tau \right)}\\ & =p\left(A_t\right)A_t+q_{\tau }\left(q_{\tau -1}\right(J_{t-1}^{(\tau -1)}\left)J_{t-1}^{(\tau -1)}\right)q_{\tau -1}\left(J_{t-1}^{(\tau -1)}\right)J_{t-1}^{(\tau -1)}\\ & =⋮\\ & =p\left(A_t\right)A_t+h(\tau ,A_{t-\tau })A_{t-\tau },\end{array}$(9) ), yields (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ), where $X_t$ represents the mature population at time $t$, $A_t$.

Remark 4.1

Consider (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ).

1. Equation (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) depends explicitly on the delay, because $\text{β}=(1+D{)}^{\tau +1}$.

2. For $rC=c,D=d$, $\tau =0$, $\text{β}=(1+D)$, (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) reduces to the classical Beverton–Holt model.

3. Equation (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) is in the form of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) with $p\left(X_t\right)$ given in (18(18) $p\left(X_t\right)=\frac{1}{1+d+c{X}_t}$(18) ), $g\left(X_t\right)=r$, and $\tilde{p}(\tau ,X_t)$ given in (22(22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) ).

4. Equation (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) satisfies conditions (H1)–(H3) with $X_{crit}<0$. Additionally, $\underset{X\to \mathrm{\infty }}{lim}h(\tau ,X)=0$.

5. $h(\tau ,X)=\tilde{p}(\tau ,X)r$ is decreasing in $\tau$.

6. The map $H$ defined in (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) is increasing in both variables, that is, (25) $\frac{\mathrm{\partial }H(u,v)}{\mathrm{\partial }u}=\frac{(1+d)}{(1+d+cu{)}^2},\frac{\mathrm{\partial }H(u,v)}{\mathrm{\partial }v}=\frac{r{D}^2\text{β}}{(D\text{β}+(\text{β}-1)rCv{)}^2}.$(25)

Remark 4.2

Note that (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) differs from the delayed Beverton–Holt model introduced in [30]: (26) $X_{t+1}=\frac{1}{1+d+c{X}_t}(X_t+\frac{rd}{\mathrm{d}\gamma +(\gamma -1)c{X}_{t-\tau +1}}{X}_{t-\tau +1}),\gamma =(1+d{)}^{\tau }.$(26) Although both models (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) and (26(26) $X_{t+1}=\frac{1}{1+d+c{X}_t}(X_t+\frac{rd}{\mathrm{d}\gamma +(\gamma -1)c{X}_{t-\tau +1}}{X}_{t-\tau +1}),\gamma =(1+d{)}^{\tau }.$(26) ) include delay solely in the growth contribution, the models are based on different assumptions regarding  the survival of the population during the delay period. More precisely, the multiplication of $s_{t+1,t}=\frac{1}{1+d+c{X}_t}$ in (26(26) $X_{t+1}=\frac{1}{1+d+c{X}_t}(X_t+\frac{rd}{\mathrm{d}\gamma +(\gamma -1)c{X}_{t-\tau +1}}{X}_{t-\tau +1}),\gamma =(1+d{)}^{\tau }.$(26) ) seems to imply that in the last time interval before a cohort of immature individuals joins the mature population, the cohort of immature individuals are exposed to competition with the already mature individuals. Furthermore, since $s_{t+1,t}$ includes the death rate $d$, one may argue that (26(26) $X_{t+1}=\frac{1}{1+d+c{X}_t}(X_t+\frac{rd}{\mathrm{d}\gamma +(\gamma -1)c{X}_{t-\tau +1}}{X}_{t-\tau +1}),\gamma =(1+d{)}^{\tau }.$(26) ) double counts the death rate $d>0$ in the last period $(t,t+1)$. Thus, (26(26) $X_{t+1}=\frac{1}{1+d+c{X}_t}(X_t+\frac{rd}{\mathrm{d}\gamma +(\gamma -1)c{X}_{t-\tau +1}}{X}_{t-\tau +1}),\gamma =(1+d{)}^{\tau }.$(26) ) may model species where immature individuals are separated from mature individuals, but have to overcome obstacles, such as travelling to the mature individuals in their last period, accounting for the additional survival factor $s_{t+1,t}$.

In contrast, in the derivation of (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ), it was assumed that each cohort of immature individuals only competes with other individuals of the same cohort until the individuals in this cohort mature. Furthermore, (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) allows for more general assumptions regarding the survival probability of immature individuals, as it allows for a different intraspecific competition rate among mature ($c$) and immature ($C$) individuals, as well as, different death rates ($d$ versus $D$). This extension is realistic, as immature individuals are often at higher risk and their survival probability is hence usually smaller.

For the remainder of this subsection, we consider (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) for $t\in {\mathbb{N}}_0$, $\tau \in \mathbb{N}=\{1,2,3,\dots ,\}$, and initial conditions ${\mathbf{X}}_{\mathbf{0}}=(X_0,X_1,\dots ,X_{\tau })$ with $X_i\ge 0$ for $i\in \{0,1,\dots ,\tau \}$.

By Proposition 3.1, solutions with non-negative initial conditions remain non-negative and, if additionally at least one initial condition is positive, then there exists $T>0$ such that $X_t>0$ for all $t>T$. Furthermore, by Proposition 3.2, solutions remain bounded.

Since $g\left(X\right)=r$ is constant, Corollary 3.8 applies so that $\overline{\tau }={\tau }_c$ with (27) ${\tau }_c=\frac{\ln \left(\frac{r(1+d)}{d(1+D)}\right)}{\ln (1+D)}.$(27) It is easy to check that (28) $\tau >{\tau }_c⟺r<(1+D{)}^{\tau +1}\frac{\mathrm{d}}{1+d}=\text{β}\frac{\mathrm{d}}{1+d}.$(28) By using Theorem 3.5 on (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ), we have the following result.

Theorem 4.3

Consider (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ). If $\tau >{\tau }_c$, then $X_0^{\ast }\equiv 0$ is the only non-negative equilibrium and $\underset{t\to \mathrm{\infty }}{lim}X\left(t\right)=X_0^{\ast }$.

Lemma 4.4

Consider (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ). If $\tau <{\tau }_c$, then there cannot be a solution $X_t$ of (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) with at least one initial condition that is positive that converges to zero.

Proof.

Since ${\tau }_c>0$, $1-p\left(0\right)<h(\tau ,0)$, and therefore, given any $ϵ>0$, there exists $\delta >0$ such that $p\left(Y\right)+h(\tau ,Y)>1+ϵ$ for all $Y\in [0,\delta ]$.

For the sake of contradiction, suppose there exists a solution with $\underset{t\to \mathrm{\infty }}{lim}{X}_t=0$. This implies that there exists $T>0$ such that $|X_{t-\tau }|<\delta$ for all $t\ge T$. Fix a $t$ so that $t>T$. Then, $X_{t-\tau },X_{t-\tau +1},\dots ,X_t<\delta$ and $X_{min}:=min\{X_{t-\tau },X_{t-\tau +1},X_{t-\tau +2},\dots ,X_t\}>0$ is well-defined, and since $F_1$ and $F_2$ are increasing, $\begin{array}{rl}{X}_{t+1}& =F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)\ge F_1\left(X_{min}\right)+F_2\left(X_{min}\right)=\left(p\right(X_{min})+h(\tau ,X_{min}\left)\right)X_{min}\\ & >(1+ϵ)X_{min}>X_{min}.\end{array}$This implies that the next-iterate will always remain above $X_{min}$, violating the assumption that the solution converges to zero.

The following theorem provides the global dynamics in the case when $\tau <\overline{\tau }={\tau }_c$.

Theorem 4.5

Consider (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ). If $\tau <{\tau }_c$, then there exists a unique positive equilibrium $X_1^{\ast }$ that is decreasing in $\tau$. Furthermore, $\underset{t\to \mathrm{\infty }}{lim}X\left(t\right)=X_1^{\ast }$ for all solutions with at least one positive initial condition.

Proof.

If $\tau <{\tau }_c$, then $h(\tau ,0)>h({\tau }_c,0)=1-p\left(0\right)$ so that, by Proposition 3.4 (c), there exists a unique positive equilibrium $X_1^{\ast }$ that is, by Lemma 3.7, decreasing in $\tau$. In fact, an explicit expression for $X_1^{\ast }$ can easily be obtained (for example using Mathematica [22]), $\begin{array}{rl}{X}_1^{\ast }& =\frac{-\left[r\mathrm{d}C\right(\text{β}-1)+cD(\text{β}-r\left)\right]}{2rcC(\text{β}-1)}\\ & +\frac{\sqrt{\left(r\mathrm{d}C\right(\text{β}-1)+cD(\text{β}-r){)}^2+4rcCD(\text{β}-1\left)\right(r(1+d)-\mathrm{d}\text{β})}}{2rcC(\text{β}-1)}.\end{array}$To show the global asymptotic stability of this equilibrium, we show that there exists $b>1$ such that $|H(X_t,X_{t-\tau })-X_1^{\ast }|<bmax\{|X_t-X_1^{\ast }|,|X_{t-\tau }-X_1^{\ast }|\}$ and then apply [20, Theorem 2]. To simplify the notation, we let $u:=X_t$ and $v:=X_{t-\tau }$.

Note that by (28(28) $\tau >{\tau }_c⟺r<(1+D{)}^{\tau +1}\frac{\mathrm{d}}{1+d}=\text{β}\frac{\mathrm{d}}{1+d}.$(28) ), $\tau <{\tau }_c$ implies that $r>\text{β}\frac{\mathrm{d}}{1+d}$. Thus, $\begin{array}{rl}H(u,v)-X_1^{\ast }& =H(u,v)-H(X_1^{\ast },X_1^{\ast })\\ & =\frac{u}{1+d+cu}+\frac{rDv}{D\text{β}+(\text{β}-1)rCv}-\frac{X_1^{\ast }}{1+d+c{X}_1^{\ast }}-\frac{rD{X}_1^{\ast }}{D\text{β}+(\text{β}-1)rC{X}_1^{\ast }}\\ & =\frac{(1+d)(u-X_1^{\ast })}{(1+d+cu)(1+d+c{X}_1^{\ast })}+\frac{r{D}^2\text{β}(v-X_1^{\ast })}{(D\text{β}+(\text{β}-1\left)rCv\right)(D\text{β}+(\text{β}-1\left)rC{X}_1^{\ast }\right)}.\end{array}$Solving $H(X_1^{\ast },X_1^{\ast })=X_1^{\ast }$ for $rD,$ it follows that $rD=\frac{(d+c{X}_1^{\ast })(D\text{β}+(\text{β}-1\left)rC{X}_1^{\ast }\right)}{(1+d+c{X}_1^{\ast })}$. Substituting this expression for $rD$ into the second term, we obtain $\begin{array}{rl}& \left|H\right(u,v)-H(X_1^{\ast },X_1^{\ast }\left)\right|\\ & =|\frac{\left(1+d\right)\left(u-X_1^{\ast }\right)}{\left(1+d+cu\right)\left(1+d+c{X}_1^{\ast }\right)}+\frac{\left(d+c{X}_1^{\ast }\right)\left(D\text{β}+\right(\text{β}-1\left)rC{X}_1^{\ast }\right)D\text{β}\left(v-X_1^{\ast }\right)}{\left(D\text{β}+\right(\text{β}-1\left)rCv\right)\left(D\text{β}+\right(\text{β}-1\left)rC{X}_1^{\ast }\right)\left(1+d+c{X}_1^{\ast }\right)}|\\ & \le \left|\frac{(1+d)(u-X_1^{\ast })}{(1+d+cu)(1+d+c{X}_1^{\ast })}\right|+\left|\frac{(d+c{X}_1^{\ast })D\text{β}(v-X_1^{\ast })}{(D\text{β}+(\text{β}-1\left)rCv\right)(1+d+c{X}_1^{\ast })}\right|\\ & \le bmax\{|u-X_1^{\ast }|,|v-X_1^{\ast }|\},\end{array}$

for $b=max\{\frac{1+d}{1+d+cϵ},\frac{D\text{β}}{D\text{β}+(\text{β}-1)rCϵ}\}$. By Lemma 4.4, the stable manifold of $X_0^{\ast }$ does not intersect the first quadrant. Therefore, by Proposition 3.3, there exists $ϵ>0$ such that $ϵ\le min\{u,v\}$ and so $b<1$. Now, the result follows by [20, Theorem 2].

Theorems 4.3 and 4.5 describe the global dynamics of (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ). Some solutions of (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) for the two different cases, $0<\tau \le {\tau }_c$ and $\tau >{\tau }_c$, are illustrated in Figure 2.

Figure 2. Solutions of (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) with parameter values $d=0.2$, $c=1$, $C=0.5$, $D=0.5$, and $r=2.5$ for different positive initial conditions. In (a), $\tau =3<{\tau }_c=5.68$. As predicted by Theorem 4.5, for $\tau <{\tau }_c$, all solutions converge to the unique positive equilibrium. In (b), $\tau =6>{\tau }_c$. Solutions converge to the only non-negative equilibrium $X_0^{\ast }=0$, as predicted by Theorem 4.3, although the convergence appears to be slow.

Remark 4.6

The global dynamics for (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) presented in Theorem 4.5 imply that the population converges to the unique positive equilibrium if the delay remains below the critical threshold in (27(27) ${\tau }_c=\frac{\ln \left(\frac{r(1+d)}{d(1+D)}\right)}{\ln (1+D)}.$(27) ) and goes extinct if the delay exceeds ${\tau }_c$. This behaviour is consistent with the global dynamics of the delayed Beverton–Holt model derived in [30]. It is also consistent with the global asymptotic stability results presented in [3] for the alternative delayed logistic differential equations model derived  there.

4.2. Ricker adult population model

We now consider the following exponential population recurrence (29) $X_{t+1}=e^{r-d-c{X}_t}{X}_t,$(29) with positive initial condition $X_0>0$, where $r>0$ determines the reproduction rate and $d>0$ is the death rate. Choosing $\tilde{r}=r-d$ and $K=\frac{r-d}{c}$ yields the classical Ricker model $X_{t+1}=X_t{e}^{\tilde{r}(1-\frac{X_t}{K})}.$Model (29(29) $X_{t+1}=e^{r-d-c{X}_t}{X}_t,$(29) ) can be expressed as (30) $X_{t+1}=e^{-d-c{X}_t}{e}^r{X}_t=e^{-d-c{X}_t}{X}_t+e^{-d-c{X}_t}(e^r-1)X_t.$(30) We can identify the survival probability from time $t$ to $t+1$ as a Poisson-type probability (31) $p\left(X_t\right)=\hat{p}\left(X_t\right)=e^{-d-c{X}_t}$(31) and the growth contribution as $g\left(X_t\right)=e^r-1=:R>0$.

Introducing a delay in the growth contribution in the sense of Section 2, we consider (32) $X_{t+1}=e^{-d-c{X}_t}{X}_t+\tilde{p}(\tau ,X_{t-\tau })R{X}_{t-\tau },$(32) with nonnegative initial conditions $X_i\ge 0$ for $i=\{0,1,\dots ,\tau \}$. To determine the survival probability of a cohort of immature individuals from time $t-\tau$ to time $t+1$, we make the same assumptions as in Section 4.1. More precisely, we assume that

1. The death rate of each cohort of immature individuals is $D>0$.

2. Each cohort of immature individuals compete only within their own cohort until the cohort reaches maturity. The competition factor is assumed to be $C>0$.

These assumptions imply that only the mature individuals have a Ricker survival probability, where as the cohorts of immature individuals have a survival probability that is determined by the Beverton–Holt recurrence.

Since the assumptions on the survival probabilities (i) and (ii) for the cohort of immature individuals are identical to the corresponding assumptions in Section 4.1, the survival probability $\tilde{p}(\tau ,X_{t-\tau })$ is given by (22(22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) ) with $r$ in (22(22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) ) replaced by $R=e^r-1$.

The delayed Ricker-adult model is then given by (33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) where $\text{β}=(1+D{)}^{\tau +1},d,c,D,R,C>0$ and non-negative initial conditions, i.e. $X_i\ge 0$ for $i\in \{0,1,\dots ,\tau \}$. Since immature individuals have a Beverton–Holt survival probability, model (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) does not simplify to the classical Ricker model for $\tau =0$. Equation (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) is in the form of (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) with $p\left(X_t\right)$ given in (31(31) $p\left(X_t\right)=\hat{p}\left(X_t\right)=e^{-d-c{X}_t}$(31) ), $g\left(X_t\right)=R=e^r-1$, and $\tilde{p}(\tau ,X_t)$ given in (22(22) $\tilde{p}(\tau ,X_{t-\tau })=\frac{D}{D(1+D{)}^{\tau +1}+((1+D{)}^{\tau +1}-1)Cr{X}_{t-\tau }}.$(22) ) with $r$ replaced by $R$. This model differs from the existing delayed Ricker population models discussed in [18, 21, 26, 32, 35].

As for the adult Beverton–Holt model (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ), the adult Ricker model (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) can be associated with an age-structured population model of the form ${\overrightarrow{Y}}_{t+1}=B_t{\overrightarrow{Y}}_t$ with $B_t$ given in (8(8) $B_t=\left[\begin{array}{cccccc}p\left(A_t\right)& 0& 0& \dots & 0& q_{\tau }\left(J_t^{\left(\tau \right)}\right)\\ q_0\left(A_t\right)& 0& 0& \dots & 0& 0\\ 0& q_1\left(J_t^{\left(1\right)}\right)& 0& \dots & 0& 0\\ ⋮& ⋮& ⋮& ⋮⋮⋮& ⋮& ⋮\\ 0& 0& 0& \dots & q_{\tau -1}\left(J_t^{(\tau -1)}\right)& 0\end{array}\right],$(8) ). Since the assumptions for the immature individuals follow the same model assumptions as in the previous example, the functions $q_i\left(w\right)$ for $i=0,1,\dots ,\tau$ are the same as in (24(24) $\begin{array}{rl}p\left(A\right)& \stackrel{\left(18\right)}{=}\frac{1}{1+d+cA},q_0\left(A\right)=\frac{r}{1+D+CrA},\\ q_i\left(w\right)& \stackrel{\left(20\right)}{=}\frac{1}{1+D+Cw},\text{for}1\le i\le \tau .\end{array}$(24) ) (replacing $r$ by $R$). The mature individuals, however, have Ricker survival probabilities, so that $p\left(A\right)=e^{-d-cA}$.

The following properties of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) are easily verified and their proofs are therefore omitted.

Remark 4.7

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ).

1. Equation (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) satisfies conditions (H1)–(H3) with $X_{crit}\le 0$ and $\underset{X\to \mathrm{\infty }}{lim}h(\tau ,X)=0$.

2. $h(\tau ,X)=\tilde{p}(\tau ,X)R$ is decreasing in $\tau$.

3. The map $H_R(u,v)$ defined in (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) is strictly increasing in $u$ for $0<cu<1$ and $v\ge 0$, and strictly decreasing for $cu>1$ and $v\ge 0$. Furthermore, $H_R(u,v)$ is strictly increasing in $v$ for all $v>0$ and $u\ge 0$. More precisely, (34) $\frac{\mathrm{\partial }{H}_R(u,v)}{\mathrm{\partial }u}=e^{-d-cu}(1-cu),\frac{\mathrm{\partial }{H}_R(u,v)}{\mathrm{\partial }v}=\frac{\text{β}{D}^{2}R}{(\text{β}D+(\text{β}-1)RCv{)}^2}.$(34)

4. By Proposition 3.1, solutions of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) with non-negative initial conditions remain non-negative for $t\ge 0$. In fact, solutions of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) with positive initial conditions remain positive for $t\ge 0$.

5. By Proposition 3.2, solutions with non-negative initial conditions remain bounded.

For the remainder of this subsection, we consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) for $t\in {\mathbb{N}}_0$, $\tau \in \mathbb{N}=\{1,2,3,\dots ,\}$, and initial conditions ${\mathbf{X}}_{\mathbf{0}}=(X_0,X_1,\dots ,X_{\tau })$ with $X_i\ge 0$ for $i\in \{0,1,\dots ,\tau \}$.

Since $g\left(X\right)=R$ is constant, Corollary 3.8 applies, so that $\overline{\tau }={\tau }_c$ with (35) ${\tau }_c=\frac{\ln \left(\frac{R}{(1+D)(1-e^{-d})}\right)}{\ln (1+D)}.$(35) It is easy to check that (36) $\tau >{\tau }_c⟺R<\text{β}(1-e^{-d}).$(36) The next result follows from Theorem 3.5.

Theorem 4.8

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau >{\tau }_c$, then the only non-negative equilibrium of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) is $X_0^{\ast }=0$ and $\underset{t\to \mathrm{\infty }}{lim}{X}_t=X_0^{\ast }$.

Theorem 4.9

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau <{\tau }_c$, then there exists a unique positive equilibrium, $X_1^{\ast }$, with value that is decreasing in $\tau$.

Proof.

If $\tau <{\tau }_c$, then $1-p\left(0\right)=h({\tau }_c,0)<h(\tau ,0)$, so that by Proposition 3.4 (c), there exists a unique positive equilibrium $X_1^{\ast }$ that is decreasing in $\tau$ by Lemma 3.7.

To address the global dynamics when $\tau <{\tau }_c$, we distinguish between two cases: (i) $c{X}_1^{\ast }<2$ and (ii) $c{X}_1^{\ast }\ge 2$ and approach the problem slightly differently in each case. We also define $\begin{array}{rl}{F}_1\left(X\right):=p\left(X\right)X& =e^{-d-cX}X,F_2\left(X\right):=h(\tau ,X)X=\frac{DRX}{m\left(X\right)},\\ & m\left(X\right):=D\text{β}+(\text{β}-1)CRX,\end{array}$so that $H_R(u,v)=F_1\left(u\right)+F_2\left(v\right)$.

Lemma 4.10

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau <{\tau }_c$, then there cannot be a solution $X_t$ of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) with at least one positive initial condition that converges to zero.

Proof.

Note that if ${\tau }_c>0$, then $1-p\left(0\right)<h(\tau ,0)$ for $0<\tau <{\tau }_c$. This proof is the same as the proof of Lemma 4.4, if one takes $0<\delta <c^{-1}$.

Theorem 4.11

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). Assume $\tau <{\tau }_c$. If $X_t$ is a solution of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) with at least one positive initial condition, then there exists $ϵ>0$ such that $\underset{t\to \mathrm{\infty }}{lim inf}{X}_t>ϵ$.

Proof.

By Lemma 4.10, the stable manifold of $X_0^{\ast }$ does not intersect the interior of the first quadrant, and so the result follows by Proposition 3.3.

The following two lemmas will be useful.

Lemma 4.12

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau <{\tau }_c$ and $c{X}_1^{\ast }<2$, then $\phi \left(X\right):=p\left(X_1^{\ast }\right)(X-X_1^{\ast })+p\left(X\right)X$ is strictly increasing for $X\ge 0$.

Proof.

For $p\left(X\right)=e^{-d-cX}$, the derivatives of $\phi \left(X\right)$ are given by $\begin{array}{rl}{\phi }^{\prime }\left(X\right)& =p\left(X_1^{\ast }\right)+p^{\prime }\left(X\right)X+p\left(X\right)=p\left(X_1^{\ast }\right)+p\left(X\right)(1-cX)\\ {\phi }^{″}\left(X\right)& =p^{″}\left(X\right)X+2p^{\prime }\left(X\right)=cp\left(X\right)(cX-2).\end{array}$Since ${\phi }^{″}\left(X\right)=0$ for $X>0$ if and only if $cX=2$, $\underset{X>0}{min}{\phi }^{\prime }\left(X\right)\ge min\left\{{\phi }^{\prime }\right(0),\underset{X\to \mathrm{\infty }}{lim}{\phi }^{\prime }(X\left)\right\}>0,$because ${\phi }^{\prime }\left(0\right)=p\left(X_1^{\ast }\right)+p\left(0\right)>0$, and since $Xp\left(X\right)\to 0$ as $X\to \mathrm{\infty }$, $\underset{X\to \mathrm{\infty }}{lim}{\phi }^{\prime }\left(X\right)>0$. Thus, ${\phi }^{\prime }\left(X\right)>0$ for all $X\ge 0$, completing the claim.

Lemma 4.13

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau <{\tau }_c$ and $c{X}_1^{\ast }<2$, then (37) $|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|\le p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$(37) where equality holds if and only if $X=X_1^{\ast }$.

Proof.

Since $p\left(X\right)$ is decreasing, $\text{sign}\left(p\right(X)-p(X_1^{\ast }\left)\right)=-\text{sign}(X-X_1^{\ast })$ so that $\begin{array}{rl}|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|& =|\left(p\right(X)-p(X_1^{\ast }\left)\right)X+p\left(X_1^{\ast }\right)(X-X_1^{\ast })|\\ & =max\left\{q_1\right(X)-q_2(X),q_2(X)-q_1(X\left)\right\},\end{array}$where $q_1\left(X\right):=|p\left(X\right)-p\left(X_1^{\ast }\right)|X\ge 0\text{and}{q}_2\left(X\right):=p\left(X_1^{\ast }\right)|X-X_1^{\ast }|\ge 0.$Clearly, if $q_2\left(X\right)-q_1\left(X\right)\ge q_1\left(X\right)-q_2\left(X\right)$, then (37(37) $|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|\le p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$(37) ) is satisfied, since $|F_1\left(X\right)-F_2\left(X\right)|=q_2\left(X\right)-q_1\left(X\right)\le q_2\left(X\right)=p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$where equality holds if and only if $X=X_1^{\ast }$.

It is left to show that if $q_1\left(X\right)-q_2\left(X\right)>q_2\left(X\right)-q_1\left(X\right)$, then (37(37) $|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|\le p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$(37) ) holds. Showing (37(37) $|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|\le p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$(37) ) holds is equivalent to showing that $q_1\left(X\right)\le 2q_2\left(X\right)$, since $p\left(X_1^{\ast }\right)|X-X_1^{\ast }|=q_2\left(X\right)$. Note that in this case, $X\ne X_1^{\ast }$, otherwise $q_1\left(X\right)=q_2\left(X\right)$ and the previous case applies.

First consider $X>X_1^{\ast }$. Then, $q_1\left(X\right)\le 2q_2\left(X\right)⟺\left(p\right(X_1^{\ast })-p(X\left)\right)X\le 2p\left(X_1^{\ast }\right)(X-X_1^{\ast }).$Rearranging terms yields the equivalent condition $p\left(X_1^{\ast }\right)(X-X_1^{\ast })+p\left(X\right)X\ge p\left(X_1^{\ast }\right)X_1^{\ast }$that is satisfied, recognizing that the left-hand side is $\phi \left(X\right)$ and the right-hand side is $\phi \left(X_1^{\ast }\right)$, for $\phi \left(X\right)$ defined in Lemma 4.12.

Next consider $X<X_1^{\ast }$. Then (37(37) $|F_1\left(X\right)-F_1\left(X_1^{\ast }\right)|\le p\left(X_1^{\ast }\right)|X-X_1^{\ast }|,$(37) ) is equivalent to $\left(p\right(X)-p(X_1^{\ast }\left)\right)X\le 2p\left(X_1^{\ast }\right)(X_1^{\ast }-X)⟺\underset{=\phi \left(X_1^{\ast }\right)}{\underbrace{p\left(X_1^{\ast }\right)X_1^{\ast }}}\ge \underset{=\phi \left(X\right)}{\underbrace{p\left(X_1^{\ast }\right)(X-X_1^{\ast })+p\left(X\right)X}},$for $\phi \left(X\right)$ given in Lemma 4.12. Since $\phi \left(X\right)$ is strictly increasing, by Lemma 4.12, the inequality holds, completing the proof.

Theorem 4.14

Consider (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). If $\tau <{\tau }_c$, and $c{X}_1^{\ast }<2$, then $X_1^{\ast }$ is globally asymptotically stable for solutions with at least one positive initial condition.

Proof.

For $X_t,X_{t-\tau }>0$, we have $\begin{array}{rl}|X_{t+1}-X_1^{\ast }|& =|F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)-F_1\left(X_1^{\ast }\right)-F_2\left(X_1^{\ast }\right)|\\ & \le |F_1\left(X_t\right)-F_1\left(X_1^{\ast }\right)|+|F_2\left(X_{t-\tau }\right)-F_2\left(X_1^{\ast }\right)|\\ & \stackrel{Lemma4.13}{\le }p\left(X_1^{\ast }\right)|X_t-X_1^{\ast }|+|F_2\left(X_{t-\tau }\right)-F_2\left(X_1^{\ast }\right)|\\ & =p\left(X_1^{\ast }\right)|X_t-X_1^{\ast }|+\frac{D^{2}R\text{β}}{m\left(X_{t-\tau }\right)m\left(X_1^{\ast }\right)}|X_1^{\ast }-X_{t-\tau }|\\ & \le \left(p\right(X_1^{\ast })+\frac{D^{2}R\text{β}}{m\left(X_{t-\tau }\right)m\left(X_1^{\ast }\right)})max\{|X_t-X_1^{\ast }|,|X_1^{\ast }-X_{t-\tau }|\}\\ & =\left(p\right(X_1^{\ast })+\frac{D\text{β}}{m\left(X_{t-\tau }\right)}(1-p\left(X_1^{\ast }\right)\left)\right)max\{|X_t-X_1^{\ast }|,|X_1^{\ast }-X_{t-\tau }|\},\\ & \le \left(p\right(X_1^{\ast })+\frac{D\text{β}}{m\left(ϵ\right)}(1-p\left(X_1^{\ast }\right)\left)\right)max\{|X_t-X_1^{\ast }|,|X_1^{\ast }-X_{t-\tau }|\},\end{array}$where we used $\frac{DR}{m\left(X_1^{\ast }\right)}=1-p\left(X_1^{\ast }\right)$. Since by Theorem 4.11, there exists $ϵ>0$ such that $X_{t-\tau }\ge ϵ>0$ for all sufficiently large $t$ and $b:=p\left(X_1^{\ast }\right)+\frac{D\text{β}}{m\left(ϵ\right)}(1-p(X_1^{\ast }\left)\right)=\left(\frac{\left(\right(\text{β}-1\left)CRϵ\right)p\left(X_1^{\ast }\right)+D\text{β}}{\left(\right(\text{β}-1\left)CRϵ\right)+D\text{β}}\right)<1$, the result follows by [20, Theorem 2].

By Theorem 4.14, when $X_1^{\ast }$ exists, it is globally asymptotically stable if $c{X}_1^{\ast }<2$. Next we discuss what we can say about the stability of $X_1^{\ast }$ when $c{X}_1^{\ast }\ge 2.$ The proof of the following theorem is provided in Appendix A1.

Theorem 4.15

Assume $\tau <{\tau }_c$ and $c{X}_1^{\ast }\ge 2$. Let $z_r$ be the largest positive root of the polynomial $N_q\left(z\right)=a_0+a_{1}z+a_2{z}^2+a_3{z}^3$ and let ${\omega }_r$ be the unique positive root of the related polynomial ${\hat{N}}_q\left(\omega \right)={\hat{a}}_0+a_{1}z+a_2{z}^2+a_3{z}^3$, where $\begin{array}{rl}{a}_3& =c(\text{β}-1{)}^2{C}^2{R}^2,\\ a_2& =(\text{β}-1)CR\left(2\text{β}\right(cD-CR)+R(2C-cD\left)\right),\\ a_1& =D(\text{β}\left[2\text{β}\right(cD-CR)+R(2C-cD\left)\right]-{\text{β}}^{2}cD-\text{β}CR(\text{β}-1)-RC(\text{β}-1\left)\right(\text{β}-R\left)\right),\\ a_0& =2\text{β}{D}^2(R-\text{β})\\ {\hat{a}}_0& =-2{\text{β}}^2{D}^2.\end{array}$

(i)	If $c{X}_1^{\ast }\in [2,c{z}_r)$, then $X_1^{\ast }$ is locally asymptotically stable.
(ii)	If $\tau$ is odd and $c{X}_1^{\ast }>max\{2,c{z}_r\}$, then $X_1^{\ast }$ is unstable.
(iii)	If $\tau$ is even and $c{X}_1^{\ast }>max\{2,c{\omega }_r\}$, then $X_1^{\ast }$ is unstable.

To determine whether the conditions of Theorem 4.15(i) or (ii) are satisfied, it suffices to check that $N_q\left(X_1^{\ast }\right)<0$ or $N_q\left(X_1^{\ast }\right)>0$, respectively and for Theorem 4.15(iii), it suffices to check that ${\hat{N}}_q\left(X_1^{\ast }\right)>0$. Figure 3 illustrates the predictions of Theorem 4.15. For the parameters used for Figure 3(a ,b), $c{X}_1^{\ast }\in [2,c{z}_r)$, andas Theorem 4.15(i) predicts, $X_1^{\ast }$ is locally asymptotically stable. Although we only showed a few representative solutions, solutions for many more randomly generated initial conditions (not shown) converge to the positive equilibrium, suggesting that $X_1^{\ast }$ is indeed globally asymptotically stable. For the parameters used for Figure 3(c) –(e), $X_1^{\ast }>z_r>2c^{-1}$ and $\tau =3$ is odd. The characteristic polynomial in (13(13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) ) in this case has a root $\lambda =-1.025$ that is outside the unit-circle, consistent with Theorem 4.15(ii), since $X_1^{\ast }$ is unstable. In fact, in this case, the solution shown converges to a periodic two cycle. This behaviour is representative as we simulated solutions for many different randomly generated initial conditions. For parameter values in Figure 3(f)–(h), Theorem 4.15 is inconclusive since $c{X}_1^{\ast }\in (max\{2,c{z}_r\},max\{2,c{\omega }_r\left\}\right)$. However, simulations of the solution suggest that $X_1^{\ast }$ is locally asymptotically stable. In fact, calculating the roots of the characteristic polynomial (13(13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) ) yields $\lambda =0.939$ and $-0.526\pm i0.843$, so that $|\lambda |<0.994<1$, confirming the local asymptotic stability of $X_1^{\ast }$. This highlights that if $\tau$ is odd, the sufficient condition in [9] is, in some sense, sharp. That is, reversing the inequality in [9] results in instability. However, for $\tau$ even, reversing the inequality in [9] does not necessarily imply instability as subfigures (f)–(h) show.

Figure 3. Solutions of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ). In all cases the initial conditions were generated randomly. (a)–(b) The parameter values are $d=0.1,c=1,C=0.05,D=0.1,R=10.$ For the given parameter values $\lceil {\tau }_c\rceil =47.$ In (a), $\tau =3$, resulting in $X_1^{\ast }\approx 3.770<z_r\approx 4.967$. In (b), $\tau =5$, resulting in $X_1^{\ast }=2.372<z_r\approx 3.584$. In both cases, solutions converge to $X_1^{\ast }$, as predicted by Theorem 4.15(i). In (c)–(e) parameter values are $d=0.01,D=1,c=1,C=0.001,R=16.2$ and $\tau =3$. Here, ${\tau }_c=9.6690$ and $X_1^{\ast }=3.3138>z_r\approx 3.287$. Theorem 4.15(ii) predicts that $X_1^{\ast }$ is unstable. The solution converges to a period two cycle. (f)–(h) The parameter values are $c=1,d=0.1,D=1,C=0.001,R=7.7$, $\tau =2$. Here, $X_1^{\ast }\approx 2.795>2c^{-1}=2$, $z_r\approx 2.275$, and ${\omega }_r\approx 15.152$. Thus, $X_1^{\ast }\in (max\{2c^{-1},z_r\},{\omega }_r)$ and Theorem 4.15 is inconclusive. However, the simulation shows convergence to $X_1^{\ast }$.

In summary, we obtained the global dynamics for the Ricker-adult delay model (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) when $\tau >{\tau }_c$. In that case, $X_0^{\ast }$ is globally asymptotically stable. Thus, as for the Beverton–Holt delay model, by Theorem 3.5, the species goes extinct if the delay exceeds the critical delay threshold, ${\tau }_c$, given in (35(35) ${\tau }_c=\frac{\ln \left(\frac{R}{(1+D)(1-e^{-d})}\right)}{\ln (1+D)}.$(35) ). For $\tau <{\tau }_c$, we proved the global asymptotic stability of $X_1^{\ast }>0$ if $c{X}_1^{\ast }<2$. Furthermore, when $c{X}_1^{\ast }\ge 2$, we provided conditions for the local asymptotic stability of $X_1^{\ast }$ and for its instability. Based on our simulations, we conjecture that local asymptotic stability of $X_1^{\ast }$ implies global asymptotic stability.

The derivation of (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) assumed a Beverton–Holt type survival for the cohorts of immature individuals. Alternatively, one may assume that the survival of the immature individuals also follows a Ricker form. For example, one may choose as the survival of one cohort with size $R{X}_{t-\tau }$ from time $t-\tau$ to $t+1$, (38) $\tilde{p}(\tau ,X_{t-\tau })={\left(e^{-D-CR{X}_{t-\tau }}\right)}^{\tau }=e^{-\tau (D+CR{X}_{t-\tau })}.$(38) Then, with everything else the same, (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ) becomes (39) $X_{t+1}=e^{-d-c{X}_t}{X}_t+R{e}^{-\tau (D+CR{X}_{t-\tau })}{X}_{t-\tau }.$(39) As for (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ), $p\left(u\right)$ satisfies (H1) and $\tilde{p}(\tau ,u)$, given in (38(38) $\tilde{p}(\tau ,X_{t-\tau })={\left(e^{-D-CR{X}_{t-\tau }}\right)}^{\tau }=e^{-\tau (D+CR{X}_{t-\tau })}.$(38) ), satisfies (H2). Furthermore, $h(\tau ,u)=\tilde{p}(\tau ,u)R$ satisfies (H3) with $X_{crit}\le 0$ and $h(\tau ,u)u$ is bounded. Thus, by Propositions 3.1–3.2, solutions remain non-negative for non-negative initial conditions, are bounded, and remain strictly positive for positive initial conditions. By Proposition 3.4, there exists at most one positive equilibrium $X_1^{\ast }$ and by Theorem 3.5, there exists a critical delay threshold ${\tau }_c$ such that for $\tau >{\tau }_c$, the population goes extinct. The critical delay threshold can be obtain by Remark 3.6, by solving $1-e^{-d}=1-p\left(0\right)=h(\tau ,0)=R{e}^{-D{\tau }_c}\text{so that}{\tau }_c=\frac{\ln \left(R\right)-\ln (1-e^{-d})}{D}.$

5. Conclusion

We presented a technique to formulate discrete delay population models by introducing a delay to represent the number of breeding cycles and therefore the number of time steps it takes a newborn individual to reach reproductive maturity, also referred to as ‘adulthood’. We also included a delay dependent survival probability to model the survival of immature individuals until these individuals reach maturity. This led to a class of delay discrete population models given by (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) that has biologically meaningful characteristics. For example, solutions remain non-negative for non-negative initial conditions. This property is not necessarily satisfied for delay recurrences that were derived using an Euler discretization scheme.

We demonstrated a relationship between our delay models and age-structured population models, that, after focussing on the adult population only, take the form of our delay models. In this case, the Leslie matrix is density dependent, and based on the model assumptions, each survival function is cohort dependent (not total population dependent). Under the assumption of invertible survival functions, as it is the case in our specific examples in Section 4, the density of each age-class can be derived from the adult population. Hence, our study of the equilibrium of our delay model also provides insight into the equilibrium age-distribution. However, to study the impact of the delay, the formulation of a delay recurrence was preferred.

Under realistic assumptions, formulated in (H1)–(H3), this class of delay discrete population models (7(7) $X_{t+1}=F_1\left(X_t\right)+F_2\left(X_{t-\tau }\right)=p\left(X_t\right)X_t+h(\tau ,X_{t-\tau })X_{t-\tau },$(7) ) has a unique critical delay threshold ${\tau }_c$. If the delay exceeds the critical delay threshold, then the population goes extinct. This is biologically justified, since the time until individuals reach their maturity prolongs the time until newborns can contribute to reproduction and increases the risk of individuals dying prior to reaching maturity. Also, if the delay is below the critical delay threshold, then the magnitude of the largest positive equilibrium as a function of the delay is decreasing. This implies a negative impact of increasing delay on the species even in the case when $0<\tau <{\tau }_c$.

Many of these properties differ from properties of other existing delay population models. For example, when introducing a delay solely in the fitness function and therefore considering the recurrence $X_{t+1}=X_{t}F\left(X_{t-\tau }\right)$, the magnitude of the equilibrium is independent of the delay and is in fact the same as for the non-delayed model $X_{t+1}=X_{t}F\left(X_t\right)$. Furthermore, the stability of the extinction equilibrium need not depend on the delay. For example, the trivial solution of the Pielou model $X_{t+1}=\frac{\alpha X_t}{1+\text{β}{X}_{t-k}}$ is unstable for all $\alpha >1$, independent of the delay, since one eigenvalue of the Jacobian is $\alpha$. It is noteworthy, however, that the positive equilibrium of the Pielou model loses stability as the delay increases, see [17],¹ but there is no transcritical equilibrium involving the trivial equilibrium as in our class of delay population models. Another model of such structure is the Ricker delay model that can exhibit stable cycles giving way to chaotic dynamics [18].

Even when delay is introduced solely in the reproduction, such as in [9], if a delay-dependent survival is not considered, the magnitude of a positive equilibrium is delay independent unless it is explicitly considered in the map. Most cases however fail to incorporate such delay, (see e.g. [16]), even though the survival of the delay period should indeed be delay dependent. More recently, in [30], a discrete delay Beverton–Holt model was introduced that was shown to have a delay-dependent positive equilibrium. Although the model formulation differed from the general formulation in this manuscript, several similarities were discussed.

With specific examples, we highlighted some interesting features. For example, even though the delay may be below the critical delay threshold, the extinction equilibrium may already be locally asymptotically stable, see (15(15) $p\left(X_t\right)=\frac{1}{1.1+.5X_t^2},\tilde{p}(\tau ,X_{t-\tau })=\frac{1}{1.1+\frac{1}{2}{X}_{t-\tau }^2+2\tau },g\left(X_{t-\tau }\right)=10X_{t-\tau }.$(15) ). Another example highlighted that delay can indeed stabilize a positive equilibrium, see Figure 1, before eventually reaching the critical delay threshold that predicts the extinction of the population.

In Section 4, we applied our derivation technique to obtain two specific models: the Beverton–Holt adult model and the Ricker adult model. For each of these classical models, we first separated the surviving (mature) individuals from the (newly) joining individuals. Following the procedure described in Section 2, we introduced a delay in the growth contribution accounting for a delay in the reproduction that may span over several breeding cycles. In both models, we assumed that the probability that individuals of a cohort reach maturity is described by a Beverton–Holt survival function. This led us to a Beverton–Holt delay model that differs from the previously introduced delay Beverton–Holt model in [30]. The difference between the delayed Beverton–Holt model in [30] and (23(23) $X_{t+1}=H(X_t,X_{t-\tau }):=\frac{1}{1+d+c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)Cr{X}_{t-\tau }}r{X}_{t-\tau }t\ge \tau ,$(23) ) is due to different model assumptions. For example, in [30], the immature individuals had the same death rate and competition rate as the mature individuals. Instead, in this manuscript, the mature individuals are exposed to a death rate $d$ and competition rate $c$ that may differ from the death rate $D$ and competition rate $C$ of immature individuals. Furthermore, in [30], the immature individuals were assumed to be exposed to mature individuals during the last time interval before maturing. Hence, the model took the form $X_{t+1}=F_1\left(X_t\right)+F_2(X_t,X_{t-\tau })$. Instead, our class of discrete delay models is based on the assumption that a cohort of immature individuals is only exposed to a survival function dependent on its cohort's density until it reaches maturity.

Despite their different model assumptions, both models exhibit similar dynamics. We discussed the global dynamics of this new Beverton–Holt adult model. We showed that if the delay exceeds the critical delay threshold ${\tau }_c$, the population goes extinct. If instead, $\tau <{\tau }_c$, then the population converges to the unique positive equilibrium with a value that is decreasing in the delay $\tau$.

For the second model that we derived, the Ricker-adult delay population model (30(30) $X_{t+1}=e^{-d-c{X}_t}{e}^r{X}_t=e^{-d-c{X}_t}{X}_t+e^{-d-c{X}_t}(e^r-1)X_t.$(30) ), we also identified the model specific critical delay threshold ${\tau }_c$ in (35(35) ${\tau }_c=\frac{\ln \left(\frac{R}{(1+D)(1-e^{-d})}\right)}{\ln (1+D)}.$(35) ). Again, if the delay exceeds this critical delay threshold, then the population goes extinct. Otherwise, if the delay remains below this critical delay threshold and the unique positive equilibrium remains bounded by twice the inverse of the competitive factor of the mature individuals (i.e. $c{X}_1^{\ast }<2$), then all solutions with positive initial conditions converge to $X_1^{\ast }$. When $c{X}^{\ast }\ge 2$, we provided conditions that guarantee the local asymptotic stability of $X_1^{\ast }$. We conjecture that in this case, $X_1^{\ast }$ is indeed globally asymptotically stable as suggested by simulations. However, an attracting two cycle is possible when the local asymptotic stability condition is strictly violated. Even though, we only discussed two specific applications of our derivation technique in Section 4, the method described in Section 2 can be applied to a variety of other discrete single species population models.

Acknowledgments

We thank the handling editor and the reviewers for their suggestions that have lead to an improvement in the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

The research of Gail S. K. Wolkowicz was partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant.

Notes

1 Based on [17], the positive equilibrium is unstable if $\frac{\alpha -1}{\alpha }>2\cos \left(\frac{k\pi }{2k+1}\right)$, where the right-hand side approaches zero as $k\to \mathrm{\infty }$.

Appendix

Appendix A1: Proof of Theorem 4.15

Consider $0<\tau <{\tau }_c$ and $c{X}_1^{\ast }\ge 2$.

(i) To show that $X_1^{\ast }$ is locally asymptotically stable, it suffices, by [9], to show that for $H_R(u,v)$ defined in (33(33) $X_{t+1}=H_R(X_t,X_{t-\tau }):=e^{-d-c{X}_t}{X}_t+\frac{D}{D\text{β}+(\text{β}-1)RC{X}_{t-\tau }}R{X}_{t-\tau },$(33) ), $|a|+|b|<1$ where (A1) $a=\frac{\mathrm{\partial }{H}_R(X_1^{\ast },X_1^{\ast })}{\mathrm{\partial }u}=e^{-d-c{X}_1^{\ast }}(1-c{X}_1^{\ast })<0,b=\frac{\mathrm{\partial }{H}_R(X_1^{\ast },X_1^{\ast })}{\mathrm{\partial }v}=\frac{\text{β}{D}^{2}R}{m^2\left(X_1^{\ast }\right)}>0.$(A1) We therefore wish to show that $\tilde{Q}\left(X_1^{\ast }\right):=|a|+|b|=e^{-d-c{X}_1^{\ast }}|1-c{X}_1^{\ast }|+\frac{\text{β}{D}^{2}R}{m^2\left(X_1^{\ast }\right)}<1,$for $m\left(x\right)=\text{β}D+(\text{β}-1)RCx$. Note that $|1-c{X}_1^{\ast }|=c{X}_1^{\ast }-1$ since $c{X}_1^{\ast }\ge 2$, and since $X_1^{\ast }$ is an equilibrium, $e^{-d-c{X}_1^{\ast }}=p\left(X_1^{\ast }\right)=1-\frac{DR}{m\left(X_1^{\ast }\right)}>0$. Thus, we aim to show that (A2) $\tilde{Q}\left(X_1^{\ast }\right)=(c{X}_1^{\ast }-1)(1-\frac{DR}{m\left(X_1^{\ast }\right)})+\frac{D^2\text{β}R}{m^2\left(X_1^{\ast }\right)}<1$(A2) or equivalently, $\tilde{Q}\left(X_1^{\ast }\right)-1<0$. Using Mathematica [22], $\tilde{Q}\left(z\right)-1=\frac{N_q\left(z\right)}{D_q\left(z\right)}=\frac{a_0+a_{1}z+a_2{z}^2+a_3{z}^3}{(\text{β}D+(\text{β}-1)RCz{)}^2},$with (A3) $\begin{array}{rl}{a}_3& =c(\text{β}-1{)}^2{C}^2{R}^2>0\\ a_2& =(\text{β}-1)CR\left(2\text{β}\right(cD-CR)+R(2C-cD\left)\right)\\ a_1& =D(\text{β}\left[2\text{β}\right(cD-CR)+R(2C-cD\left)\right]-{\text{β}}^{2}cD-\text{β}CR(\text{β}-1)-RC(\text{β}-1\left)\right(\text{β}-R\left)\right)\\ a_0& =2\text{β}{D}^2(R-\text{β}).\end{array}$(A3) We note that (A4) $N_q\left(2c^{-1}\right)=a_0+2a_1{c}^{-1}+4a_2{c}^{-2}+8a_3{c}^{-3}=\frac{-2(\text{β}-1)CD{R}^2}{c}<0.$(A4) Since the leading coefficient of z in $N_q\left(z\right)$ is positive, this implies that there exists at least one positive root $z_r$ of $N_q\left(z\right)$ such that $z_r>2c^{-1}$.

Case 1: Assume that $R>\text{β}$. Then $a_0>0$, so that by Descartes' rule of signs, there exist either no positive root or two positive roots. By (A4(A4) $N_q\left(2c^{-1}\right)=a_0+2a_1{c}^{-1}+4a_2{c}^{-2}+8a_3{c}^{-3}=\frac{-2(\text{β}-1)CD{R}^2}{c}<0.$(A4) ), there exists two positive roots $z_{\ell }$ and $z_r$ such that $0<z_{\ell }<2c^{-1}<z_r$. Hence, for $2c^{-1}\le X_1^{\ast }<z_r$, (A4(A4) $N_q\left(2c^{-1}\right)=a_0+2a_1{c}^{-1}+4a_2{c}^{-2}+8a_3{c}^{-3}=\frac{-2(\text{β}-1)CD{R}^2}{c}<0.$(A4) ) holds and therefore, (A2(A2) $\tilde{Q}\left(X_1^{\ast }\right)=(c{X}_1^{\ast }-1)(1-\frac{DR}{m\left(X_1^{\ast }\right)})+\frac{D^2\text{β}R}{m^2\left(X_1^{\ast }\right)}<1$(A2) ) also holds.

Case 2: Assume that $R<\text{β}$. Then $a_0<0$. By Descartes' rule of signs, there exists either one or three positive roots. Furthermore, if $a_1>0$, then $a_2>0$, because $\frac{a_1}{D}=\frac{a_2\text{β}}{(\text{β}-1)CR}-({\text{β}}^{2}cD+\text{β}CR(\text{β}-1)+RC(\text{β}-1\left)\right(\text{β}-R\left)\right)<\frac{a_2\text{β}}{(\text{β}-1)CR}.$Thus, by Descartes' rule of signs, there exists exactly one positive root $z_r>2c^{-1}$, so that $\tilde{Q}\left(z\right)-1<0$ for $0\le z<z_r$. Hence, for $2c^{-1}\le X_1^{\ast }<z_r$, the desired inequality holds.

Case 3: Assume that $R=\text{β}$. If $a_1>0$, then by the argument in Case 2, $a_2>0$. This would imply, by Descartes' rule of signs, that there exists no positive real root. This however violates (A4(A4) $N_q\left(2c^{-1}\right)=a_0+2a_1{c}^{-1}+4a_2{c}^{-2}+8a_3{c}^{-3}=\frac{-2(\text{β}-1)CD{R}^2}{c}<0.$(A4) ). Thus, $a_1<0$, so that by Descartes' rule of signs, there exists exactly one positive real root $z_r$ such that the desired inequality holds for $0<z<z_r$. Hence, for $2c^{-1}\le X_1^{\ast }<z_r$, the desired inequality holds.

In all cases, as long as $2c^{-1}\le X_1^{\ast }<z_r$, where $z_r$ is the largest positive real root of the polynomial $N_q\left(z\right)=a_0+a_{1}z+a_2{z}^2+a_3{z}^3,$the desired inequality holds.

To prove the instability results, we consider the cases when $\tau$ is odd and when $\tau$ is even, separately.

(ii) First assume that $\tau$ is odd and assume that $X_1^{\ast }>max\{2c^{-1},z_r\}$. Since $X_1^{\ast }>z_r$, by the above argument, (A5) $\tilde{Q}\left(X_1^{\ast }\right)>1⟺|a|+|b|>1⟺-a+b>1,$(A5) for a and b defined in (A1(A1) $a=\frac{\mathrm{\partial }{H}_R(X_1^{\ast },X_1^{\ast })}{\mathrm{\partial }u}=e^{-d-c{X}_1^{\ast }}(1-c{X}_1^{\ast })<0,b=\frac{\mathrm{\partial }{H}_R(X_1^{\ast },X_1^{\ast })}{\mathrm{\partial }v}=\frac{\text{β}{D}^{2}R}{m^2\left(X_1^{\ast }\right)}>0.$(A1) ). Recall that, by (13(13) ${\lambda }^{\tau +1}-a{\lambda }^{\tau }-b=0,$(13) ), the characteristic polynomial is $P\left(\lambda \right)={\lambda }^{\tau +1}-a{\lambda }^{\tau }-b.$Since $P(-1)=(-1{)}^{\tau +1}-a(-1{)}^{\tau }-b=1+a-b=-(-1-a+b)\stackrel{\left(A5\right)}{<}0,$and $\underset{\lambda \to -\mathrm{\infty }}{lim}P\left(\lambda \right)=\mathrm{\infty }$, there exists a root of the characteristic polynomial, say $\overline{\lambda }$, such that $\overline{\lambda }<-1$ and therefore outside the unit-circle. Hence, $X_1^{\ast }$ is unstable.

(iii) Now assume that $\tau$ is even, $c{X}_1^{\ast }\ge 2$, and $X_1^{\ast }>{\omega }_r$.

First, we show that if $X_1^{\ast }>max\{2c^{-1},w_r\}$, then $-a-(1+b)>0$. Using Mathematica [22], we express $-a-(1+b)=\hat{Q}\left(\omega \right)=\frac{{\hat{N}}_q\left(\omega \right)}{{\hat{D}}_q\left(\omega \right)}=\frac{b_0+b_1\omega +b_2{\omega }^2+b_3{\omega }^3}{(\text{β}D+(\text{β}-1)CR\omega {)}^2},$where $b_i=a_i,i=1,2,3$ defined in (A3(A3) $\begin{array}{rl}{a}_3& =c(\text{β}-1{)}^2{C}^2{R}^2>0\\ a_2& =(\text{β}-1)CR\left(2\text{β}\right(cD-CR)+R(2C-cD\left)\right)\\ a_1& =D(\text{β}\left[2\text{β}\right(cD-CR)+R(2C-cD\left)\right]-{\text{β}}^{2}cD-\text{β}CR(\text{β}-1)-RC(\text{β}-1\left)\right(\text{β}-R\left)\right)\\ a_0& =2\text{β}{D}^2(R-\text{β}).\end{array}$(A3) ) and $b_0=-2{\text{β}}^2{D}^2<0$.

Since the leading coefficient is positive and ${\hat{N}}_q\left(\omega \right)\left(2c^{-1}\right)=\frac{-\left(2DR\right(\text{β}cD+(\text{β}-1)CR)}{c}<0,$there exists at least one root ${\omega }_r>2c^{-1}$.

By the above argument, we know that if $a_1>0$, then $a_2>0$. Hence, if $b_1>0$, then $b_2>0$. Thus, by Descartes' rule of signs, there exists exactly one positive real root ${\omega }_r>0$. Hence, if $c{X}_1^{\ast }>c{\omega }_r$, then $-a-(1+b)>0$. In this case, $P(-1)=(-1{)}^{\tau +1}-a(-1{)}^{\tau }-b=-1-a-b>0.$However, since $\underset{\lambda \to -\mathrm{\infty }}{lim}P\left(\lambda \right)=-\mathrm{\infty }$, there exists a root of the characteristic polynomial $P\left(\lambda \right)$, say $\hat{\lambda }$, such that $\hat{\lambda }<-1$, outside of the unit-circle. Thus, $X_1^{\ast }$ is unstable.