We establish a necessary and sufficient condition for global stability of the nonlinear discrete red blood cells survival model and demonstrate that local asymptotic stability implies global stability. Oscillation and solution bounds are investigated. We also show that, for different values of the parameters, the solution exhibits some time-varying dynamics, that is, if the system is moved in a direction away from stability (by increasing the parameters), then it undergoes a series of bifurcations that leads to increasingly long periodic cycles and finally to deterministic chaos. We also study the chaotic behavior of the model with a constant positive perturbation and prove that, for large enough values of one of the parameters, the perturbed system is again stable.

1. Introduction

In 1950–1980 various biological and medical phenomena were for the first time described using differential equations, among them the blood cell production [1–3].

The general approach was to describe the cell production dynamics as ğ‘î…ž=𝑝(𝑡,𝑁(𝑡))−𝑑(𝑡,𝑁(𝑡)),(1.1) where 𝑁 is the number of cells, 𝑝(𝑡,𝑁) is the cell production rate at time 𝑡 under the condition that the population size is 𝑁, and 𝑑(𝑡,𝑁) is the mortality rate. It was a common belief that the mortality is proportional to the amount of blood cells in the circulation 𝑑(𝑡,𝑁(𝑡))=𝜇𝑁(𝑡), where 𝜇 can be time dependent or not, and that higher per capita production rates correspond to smaller numbers of existing cells 𝑁. However, if there exists a constant equilibrium 𝑁, 𝑝 and 𝑑 are time independent and 𝑝(𝑁)=𝑑(𝑁), while 𝑝(𝑁)>𝑑(𝑁) for 0<𝑁<𝑁, 𝑝(𝑁)<𝑑(𝑁) for 𝑁>𝑁, then the dynamics of (1.1) is very simple: all solutions with 𝑁(0)>0 are positive and monotonically converge to the constant equilibrium. Real data suggest that such models poorly describe the oscillatory and chaotic behavior which frequently occurs in observations. As the next stage, it was suggested to introduce delay in the production term (the maturation delay): in fact, cell production in the bone marrow takes some time before the cell is released into circulation. The resulting delay equations can demonstrate oscillatory and chaotic behavior; see, for example, [1, 2, 4].

The delay differential equation 𝑁′(𝑡)=−𝜇𝑁(𝑡)+ğ‘ƒğ‘’âˆ’ğ‘žğ‘(𝑡−𝜏),𝑡≥0,(1.2) was introduced by Ważewska-Czyżewska and Lasota [3] to describe the survival of red blood cells in an animal. Here, as above, 𝑁(𝑡) is the number of the mature red blood cells at time 𝑡, 𝜇>0 is the per capita death rate (the probability of death for blood cells which currently circulate), and 𝑃>0 and ğ‘ž>0 define the red blood cell production function: 𝑃 can be described as the limit production when the number of cells tends to zero, the decay of cell production for large cell number becomes faster with the growth of ğ‘ž, and delay 𝜏 is the time required to produce a red blood cell. Equation (1.2) is considered with a nonnegative initial function and a positive initial value which describe the number of cells in the past: 𝑁(𝑡)=𝜑(𝑡),𝜑(𝑡)≥0,−𝜏≤𝑡<0,𝑁(0)=𝑁0>0.(1.3) Evidently (1.2) with initial condition (1.3) has a unique solution which is positive for any 𝑡≥0. Equation (1.2) has a positive equilibrium 𝑁 which satisfies the equation 𝑃𝑁=ğœ‡ğ‘’âˆ’ğ‘žğ‘.(1.4)

After (1.2) was introduced in [3], this equation and its nonautonomous modification 𝑁′(𝑡)=−𝜇(𝑡)𝑁(𝑡)+𝑃(𝑡)ğ‘’âˆ’ğ‘žğ‘(𝑡−𝜏(𝑡)),𝑡≥0,(1.5) where 𝜇(𝑡)>0, 𝑃(𝑡)>0, ğ‘ž>0, 𝜏(𝑡)≥0, were intensively studied [4–8], and some results were presented in the monograph [9]. Later on, several generalizations of (1.2) were investigated, for example, equations with impulses [10–13] and a distributed delay [14].

Discrete modifications of (1.2), (1.5) are much less studied; see, for example, [15]. Discretization of (1.5) can be obtained by the assumption that 𝑡−𝜏(𝑡) is a piecewise constant function, for example, 𝑡−𝜏(𝑡)=ℎ[𝑡/ℎ]: 𝑁′(𝑡)=−𝜇(𝑡)𝑁(𝑡)+𝑃(𝑡)ğ‘’âˆ’ğ‘žğ‘(ℎ[𝑡/ℎ]),(1.6) where [𝑥] is the integer part of 𝑥, and successive integration of this equation over the interval of length ℎ. After integrating (1.6) over [ğ‘›â„Ž,𝑡], where ğ‘›â„Žâ‰¤ğ‘¡â‰¤(𝑛+1)ℎ, we obtain at 𝑡≥𝑛𝑡𝑁(𝑡)=𝑁(𝑛𝑡)ğ‘’âˆ’âˆ«ğ‘¡ğ‘›â„Žğœ‡(𝑠)𝑑𝑠+î‚¸î€œğ‘¡ğ‘›â„Žğ‘ƒ(𝜏)𝑒−∫𝑡𝜏𝜇(𝑠)ğ‘‘ğ‘ î‚¹ğ‘’ğ‘‘ğœâˆ’ğ‘žğ‘(ğ‘›â„Ž),(1.7)

which after substituting 𝑡=(𝑛+1)ℎ and denoting 𝑥(𝑛)=𝑁(ğ‘›â„Ž),𝑛∈ℕ,(1.8)𝛼(𝑛)=1−𝑒−∫(𝑛+1)â„Žğ‘›â„Žğœ‡(𝑠)𝑑𝑠>0,𝑝(𝑛)=(𝑛+1)â„Žğ‘›â„Žğ‘ƒ(𝜏)𝑒−∫𝜏(𝑛+1)â„Žğœ‡(𝑠)𝑑𝑠𝑑𝜏(1.9)

can be rewritten as 𝑥(𝑛+1)=(1−𝛼(𝑛))𝑥(𝑛)+𝑝(𝑛)ğ‘’âˆ’ğ‘žğ‘¥(𝑛).(1.10) In the case when the delay is a multiple of the time step 𝑡𝜏(𝑡)=ℎℎ+ğ‘˜â„Ž,(1.11)

where {𝑥}=𝑥−[𝑥] is the fractional part of the number, the same integration process as above leads to the delay difference equation 𝑥(𝑛+1)=(1−𝛼(𝑛))𝑥(𝑛)+𝑝(𝑛)ğ‘’âˆ’ğ‘žğ‘¥(𝑛−𝑘).(1.12)

Compared to ordinary and delay differential equations, difference equations have certain advantages, and their study can be motivated by the following reasons.(1)On the one hand, compared to ordinary differential equations, difference models have much richer and more complex dynamics, including oscillation, cycles, and transition to chaos. On the other hand, due to a simpler form, they allow a relatively simple analysis of stability, bifurcations, and some other dynamics-related issues. Comparing two models, nondelay difference equation (1.10) and delay equation (1.12), we notice that already (1.10) inherits some main properties of the dynamics of original delay differential equation (1.2): existence of unstable positive solutions for some values of the parameters, sustainable oscillations, and chaotic behavior. (2)In practical implementations, discrete models are more convenient for both numerical analysis and parameter estimation from the scarce experimental data available. (3)Difference equations can be considered as a result of discretization or numerical approximation of differential equations, and their dynamics is interesting from the point of view of the properties of relevant numerical schemes.

As mentioned above, there have already been some publications on the discrete Lasota-Wazewska model. However, the present paper is either different or completes the existing results from the following points of view.(1)The model considered in the present paper is relatively simple (autonomous, without delay). However, for this equation we for the first time obtain strict solution bounds and sharp stability results, not just sufficient attractivity conditions. (2)The problems studied earlier for difference equations of this type were usually restricted to bounds, oscillation, and stability, and the analysis mainly referred to the range of parameters lower than those defining the first period-doubling bifurcation when the stable two-cycle appears. The present paper studies the whole range of parameters, including the chaotic areas. (3)For the first time, constant perturbation of the discrete Lasota-Wazewska model was introduced and the resulting period-halving bifurcations were investigated.

The paper is organized as follows. Section 2 presents known facts applied to (2.1) and contains some auxiliary results which will be applied later. In Section 3 we obtain lower and upper bounds for solutions of (2.1), present a sharp nonoscillation condition, and demonstrate that the local asymptotic stability of the positive equilibrium of (2.1) is equivalent to its global asymptotic stability. In Section 4 we investigate bifurcations of the map described by (2.1) in two parameters: 𝑝 and ğ‘ž. It is demonstrated that an addition of a small positive perturbation to the right hand side of (2.1) leads to stable behavior for ğ‘ž large enough but does not change the character of the bifurcation diagram in 𝑝. Finally, in Section 5 the results of the present paper are summarized and several open problems are outlined.

2. Preliminaries

In the present paper the difference equation 𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+ğ‘ğ‘’âˆ’ğ‘žğ‘¥(𝑛),𝑛≥0,(2.1) is considered, where 𝑝,ğ‘žâˆˆ(0,∞),𝛼∈(0,1).(2.2) By the biological interpretation (𝑥(𝑛) is the red blood cells count), we refer to positive solutions of (2.1) only; in particular, we assume that the initial condition is positive: 𝑥(0)>0.(2.3)

It is clear that the initial value problem (2.1) and (2.3) has a unique positive solution {𝑥(𝑛)}âˆžğ‘›=1 which is bounded and persistent (see Kubiaczyk and Saker [16] and Li and Cheng [17]). Since the function 𝑔(𝑥)=ğ›¼ğ‘¥âˆ’ğ‘ğ‘’âˆ’ğ‘žğ‘¥(2.4) has exactly one root on (0,∞), (2.1) has a unique positive fixed point 𝑥.

The following stability definitions follow the monograph [18].

Definition 2.1. The point 𝑥∗ is a critical point of the difference equation 𝑥𝑛+1𝑥=𝑓𝑛(2.5) with a continuously differentiable function 𝑓(𝑥) if 𝑓′(𝑥∗)=0. 𝑥 is a fixed point or an equilibrium point if 𝑓(𝑥)=𝑥. The equilibrium point is locally stable if for any 𝜀>0 there exists 𝛿>0 such that |𝑥0−𝑥|<𝛿 implies |𝑥𝑛+1−𝑥|=|𝑓𝑛(𝑥0)−𝑥|<𝜀. The fixed point 𝑥 is globally asymptotically stable if it is locally stable and limğ‘›â†’âˆžğ‘¥ğ‘›+1=limğ‘›â†’âˆžğ‘“ğ‘›(𝑥0)=𝑥 for any 𝑥0 in the domain including 𝑥; everywhere in this paper we assume that 𝑥0∈(0,∞).
Our purpose is to establish a new global stability condition for (2.1) which is both necessary and sufficient and prove that local stability implies global stability.
We will also use the following definition.

Definition 2.2. The difference equation 𝑥𝑛+1𝑥=𝑓𝑛(2.6) is persistent if there exists 𝑚>0 such that any solution with 𝑥0>0 satisfies limsupğ‘›â†’âˆžğ‘¥ğ‘›â‰¥ğ‘š.(2.7)

In the recent decades the problems associated with understanding the dynamical behavior of nonlinear difference equations have been receiving intensive attention. The asymptotic behavior of solutions of (2.1) has been proposed as a research project by Kocić and Ladas [19, Project 4.6.1]. Since then several authors have expounded on various aspects of this model. For contributions, we refer the reader to the papers by Karakostas et al. [20], Zheng et al. [21], Ma and Yu [22], Meng and Yan [23], Li and Cheng [17], Ivanov [24], Györi and Trofimchuk [25], El-Morshedy and Liz [26], Wang and Li [15], and Saker [27].

The stability of fixed points is one of the most important issues in the study of population dynamics models. In fact, it is much easier to perform local stability analysis of a fixed point than its global stability analysis. For local stability of (2.1) we see that, if ||ğ‘“î…žî€·ğ‘¥î€¸||=||1âˆ’ğ›¼âˆ’ğ‘ğ‘žğ‘’âˆ’ğ‘žğ‘¥||<1,(2.8) where 𝑓(𝑥)∶=(1−𝛼)𝑥+ğ‘ğ‘’âˆ’ğ‘žğ‘¥,(2.9) then the fixed point 𝑥 is locally stable. As a special case of the results established in the above mentioned papers for global stability, we have the following global stability conditions for the unique fixed point 𝑥 of (2.1).

We note that 𝛽 in (C7) is the unique solution of the equation ℎ2(𝑥)=𝑥 in (0,𝑥), where ℎ(𝑥)=(1−𝛼)−1ğ‘ğ‘’âˆ’ğ‘žğ‘¥.

One can see that the global stability conditions (C1)–(C9) are different from the local stability condition (2.8). The natural question now is the following: can we obtain a new global stability condition for (2.1) which improves the above conditions? Our aim in this paper is to give an affirmative answer to this question and prove that the local stability condition (2.8) is also a global stability condition for (2.1). This means that the local stability implies the global stability for (2.1) Table 1.

In general, local and global stability conditions for difference equations do not coincide and different methods should be used to investigate the conditions for global stability. On the one hand, it is possible for a model to be locally stable and fail to be globally stable when the peak (maximum) is too high and when the curve falls off too quickly. On the other hand, it is possible for a model to be globally stable and fail to be locally stable if the nonlinear function is not continuous. In fact, there is no obvious connection between local and global stability conditions, and there is still a gap between the most local and global stability tests.

Cull [28, 29] considered this problem and proved that local stability implies global stability for some discrete population models. The results in [28, 29] are based on the enveloping principle, where it is proved that the enveloping by a linear fractional function is sufficient for global stability. However, this proof is not easy for complicated models and the parameter in the enveloping function must be adjusted for each particular population model. Moreover, the map should be a unimodal function in the sense of [28], which is not applicable to the model (2.1) that we will study in this paper.

Further we will apply the following result on the equivalence of local and global stability; its exact statement is taken from [25, Proposition  7], see also [30, 31].

Lemma 2.3. Let the function 𝑔∶[ğ‘Ž,𝑏]→[ğ‘Ž,𝑏], where 𝑔∈𝐶3[ğ‘Ž,𝑏], have at most one critical point 𝑥∗ in [ğ‘Ž,𝑏]. If a unique fixed point 𝑥∈[ğ‘Ž,𝑏] of 𝑔 is locally asymptotically stable and the Schwarzian derivative 𝑔(𝑆𝑔)(𝑥)=(𝑥)ğ‘”î…ž(−3𝑥)2î‚µğ‘”î…žî…ž(𝑥)ğ‘”î…žî‚¶(𝑥)2(2.10) is negative for 𝑥∈[ğ‘Ž,𝑏], 𝑥≠𝑥∗, then 𝑥 is globally asymptotically stable.

3. Global Attractivity

In this section, we study qualitative properties of (2.1) and we establish a new global stability condition by applying Lemma 2.3. Before we state stability results, let us discuss the properties of the function 𝑓(𝑥)∶=(1−𝛼)𝑥+ğ‘ğ‘’âˆ’ğ‘žğ‘¥.(3.1) We have 𝑓(0)=𝑝>0. If ğ‘ğ‘žâ‰¤1−𝛼, then 𝑓(𝑥) is increasing for any 𝑥>0. If ğ‘ğ‘ž>1−𝛼, then 𝑓(𝑥) is decreasing for 0<𝑥<𝑥∗ and 𝑓(𝑥) is increasing for 𝑥>𝑥∗, where 𝑥∗=1ğ‘žî‚€lnğ‘ğ‘žî‚1−𝛼.(3.2) This means that (2.1) is not a population model in the sense of Cull’s definition [28], which describes positive unimodal functions increasing for 𝑥<𝑥∗ and decreasing for 𝑥>𝑥∗. Thus the results obtained, for example, in [28, 29] cannot be applied to (2.1).

Theorem 3.1. If ğ‘ğ‘žâ‰¤1−𝛼 or î‚€ğ‘ğ‘ž>1−𝛼,lnğ‘ğ‘žî‚â‰¤1−𝛼1−𝛼𝛼,(3.3) then all positive solutions of (2.1) either exceed the equilibrium 𝑥 for 𝑛≥2 or are less than 𝑥, or coincide with 𝑥 for any 𝑛≥0. If ğ‘ğ‘žâ‰¤1−𝛼 or (3.3) is satisfied, then all solutions converge monotonically to 𝑥.

Proof. Let ğ‘ğ‘žâ‰¤1−𝛼. Then 𝑓(𝑥) is increasing for any 𝑥>0, and all solutions converge monotonically to 𝑥 (thus, {𝑥(𝑛)} is increasing if 𝑥(0)<𝑥, since 𝑥(𝑛)<𝑥 implies 𝑥(𝑛)<𝑥(𝑛+1)<𝑥, and is decreasing for 𝑥(0)>𝑥).
Further, let ğ‘ğ‘ž>1−𝛼. Since 𝑓(𝑥) is decreasing for 0<𝑥<𝑥∗ and increasing for 𝑥>𝑥∗, then all positive solutions are oscillatory about the positive equilibrium whenever 𝑥<𝑥∗. In fact, if 𝑥<𝑥∗, then 𝑥(𝑛)<𝑥 implies 𝑥(𝑛+1)>𝑥, since 𝑓(𝑥)>𝑥 for any 𝑥∈[0,𝑥). There is also 𝐾>𝑥 such that 𝑓(𝐾)=𝑥 (e.g., in Figure 2(b), 𝐾≈3.9, where 𝑓(𝐾)=𝑥≈1.96). If 𝑥(𝑛)≥𝐾 then there exists 𝑗∈ℕ such that either 𝑥(𝑛+𝑗)<𝐾, which implies 𝑥(𝑛+𝑗+1)<𝑥, or all 𝑥(𝑛+𝑗)=𝑥 for 𝑗 large enough.
All solutions are nonoscillatory for 𝑥≥𝑥∗: they either exceed 𝑥 beginning with the second iteration if at least one of the inequalities 𝑥(0)>𝑥 and 𝑓(𝑥(0))>𝑥 is satisfied or are less than 𝑥 for any 𝑛, otherwise. The fact that 𝑥≥𝑥∗ is equivalent to 𝑓(𝑥∗)≥𝑥∗, or 1âˆ’ğ›¼ğ‘žî‚€lnğ‘ğ‘žî‚+1−𝛼1âˆ’ğ›¼ğ‘žâ‰¥1ğ‘žî‚€lnğ‘ğ‘žî‚î‚€1−𝛼⟺lnğ‘ğ‘žî‚â‰¤1−𝛼1âˆ’ğ›¼ğ‘ž,(3.4)completes the proof.

We note that in the nonoscillatory case we have monotone convergence to 𝑥; such monotonically convergent solutions for 𝛼=0.5,𝑝=0.6, and ğ‘ž=2 are illustrated in Figure 1.

Remark 3.2. Theorem 3.1 presents nonoscillation conditions for (2.1). Moreover, the proof yields that these conditions are sharp: if they are not satisfied, all positive solutions (which are not identically equal to 𝑥 for 𝑛 large enough) are oscillatory about 𝑥.

Remark 3.3. In the proof of Theorem 3.1 we have applied the fact that will be used in future: if 𝑓(𝑥)>0 for 𝑥>0, the map 𝑓 has the only positive equilibrium 𝑥 and is increasing on (𝑥∗,∞), where 𝑥∗<𝑥, then all trajectories with 𝑥(0)>0 are eventually monotone and converge to 𝑥.

Theorem 3.4. For any positive solution 𝑥(𝑛) of (2.1) one has liminfğ‘›â†’âˆžğ‘¥(𝑛)=limsupğ‘›â†’âˆžğ‘¥(𝑛)=𝑥(3.5) if either ğ‘ğ‘žâ‰¤1−𝛼 or (3.3) is satisfied. If none of the above holds, then liminfğ‘›â†’âˆžğ‘¥î€·ğ‘¥(𝑛)≥𝑓∗,limsupğ‘›â†’âˆžğ‘¥î€·ğ‘“î€·ğ‘¥(𝑛)≤𝑓∗,(3.6) where 𝑓(𝑥) and 𝑥 are defined in (3.1) and (3.2), respectively.

Proof. If either ğ‘ğ‘žâ‰¤1−𝛼 or (3.3) is satisfied, then by Theorem 3.1 we have monotone convergence, so we will just consider the case when ğ‘ğ‘ž>1−𝛼 and (3.3) does not hold. Since 𝑥∗ is the absolute minimum point for 𝑓(𝑥), 𝑥≥0, then all 𝑥(𝑛) but probably 𝑥(0) satisfy 𝑥(𝑛)≥𝑓(𝑥∗). We have 𝑥∗>𝑥, so the sequence 𝑥(𝑛) with 𝑥(𝑛)>𝑥∗ is decreasing until 𝑓(𝑥(𝑛))<𝑥∗. We note that, since 𝑓(𝑥) is decreasing for 𝑥∈(0,𝑥∗), then max[𝑥∈𝑓(𝑥∗),𝑥∗]𝑓𝑓𝑥(𝑥)=𝑓∗,(3.7) which completes the proof.

Remark 3.5. Let us note that in the case when we have a stable 2𝑘-cycle, sharper bounds than (3.6) can be obtained. Moreover, if ğ‘ž grows, then even for the deterministic chaos the attraction domain of the solution becomes narrow (see Figure 5(a)), for the bifurcation diagram).

Theorem 3.6. If ||1âˆ’ğ›¼âˆ’ğ‘ğ‘žğ‘’âˆ’ğ‘žğ‘¥||<1,(3.8) then the positive equilibrium point 𝑥 of (2.1) is globally stable, that is, the local stability of 𝑥 implies its global stability.

Proof. We can omit the monotonic convergence case 𝑥∗<𝑥 and assume that 𝑥∈(0,𝑥∗). Since 𝑓(0)=𝑝 and 𝑓(𝑥) is decreasing in (0,𝑥∗), then for any 𝑏≥max{𝑝,𝑥∗} the function 𝑓 maps the segment [0,𝑏] onto itself. The Schwarzian derivative of function 𝑓 defined by (3.1) is negative, (𝑆𝑓)(𝑥)=âˆ’ğ‘ğ‘ž3î€ºğ‘ğ‘ž+2ğ‘’ğ‘žğ‘¥î€»(1−𝛼)2[ğ‘’ğ‘žğ‘¥](1−𝛼)âˆ’ğ‘ğ‘ž2<0,(3.9) if 𝑥≠𝑥∗=ln(ğ‘ğ‘ž/(1−𝛼))/ğ‘ž, so by Lemma 2.3 a locally asymptotically stable equilibrium 𝑥 is also globally asymptotically stable.

To illustrate the main result in Theorem 3.6, we consider the following examples.

Example 3.7. Consider the model 1𝑥(𝑛+1)=2𝑥(𝑛)+10𝑒−2𝑥(𝑛).(3.10) Here 𝑝=10, ğ‘ž=2, and 𝛼=0.5. The fixed point of (3.10) is given by 𝑥≈1.3484. One can see that conditions (C1)–(C9) are not satisfied, so these results cannot be applied to establish stability of (3.10). Since ||1âˆ’ğ›¼âˆ’ğ‘žğ‘ğ‘’âˆ’ğ‘žğ‘¥||≈|||11−2−20𝑒−2(1.348)|||≈0.85<1,(3.11)then by Theorem 3.6 the equilibrium 𝑥 is globally stable. For illustration see Figure 2(a).

Example 3.8. Consider the model 1𝑥(𝑛+1)=2𝑥(𝑛)+50𝑒−2𝑥(𝑛).(3.12)Here 𝑝=50, ğ‘ž=2, 𝛼=1/2. We note that for the positive fixed point 𝑥≈1.965 of (3.12) ||1âˆ’ğ›¼âˆ’ğ‘žğ‘ğ‘’âˆ’ğ‘žğ‘¥||≈|||11−2−100𝑒(−2)1.965|||≈1.464>1,(3.13) thus 𝑥 is neither locally nor globally stable. Instead of a stable equilibrium, there is a stable two-cycle ≈ 2.86, 1.59, see Figure 2(b). The bifurcation diagram in Figure 4 demonstrates that 𝑝=50 is in the range of parameter 𝑝 where (2.1) has a stable 2-cycle.

Remark 3.9. We also illustrate Theorem 3.6 by applying the Coppel result [32] which claims that the single positive equilibrium 𝑥 of 𝑓(𝑥) is globally stable whenever the equation 𝑥=𝑓(𝑓(𝑥)) has the only positive solution 𝑥. For (3.10), by Theorem 3.6 the positive equilibrium point is globally stable. This is illustrated by the graph in Figure 3(a), where the curve of 𝑦=𝑓(𝑓(𝑥)) intersects 𝑦=𝑥 only once. By Theorem 3.6 for (3.12) the positive equilibrium point is unstable. Figure 3(b), illustrates that there are more than one intersection points (really, three: the middle one corresponds to the equilibrium, while the left and the right points form a stable two-cycle) of 𝑦=𝑓(𝑓(𝑥)) and 𝑦=𝑥, so 𝑥 is unstable.

4. Bifurcations and Chaos

Stability properties of (2.1) can also be illustrated using the bifurcation diagram. The critical point 𝑝𝑐 at which the unique positive equilibrium 𝑥 loses its stability (ğ‘“î…ž(𝑥)=−1) can be easily computed as 𝑝𝑐=ğ›¼î€·ğ‘žğ‘¥exp𝑥,where𝑥=2âˆ’ğ›¼ğ‘žğ›¼.(4.1) For 𝑝<𝑝𝑐, all positive solutions converge to the equilibrium point 𝑥 (see Figure 4(a)). As 𝑝 increases beyond 𝑝𝑐, there is a series of period-doubling bifurcations leading to a deterministic chaos. This means that for parameters close to the stable region, this will be a stable two-cycle and if the system is moved in a direction away from stability, by increasing the parameters, then the dynamics becomes more complex and the system undergoes a series of bifurcations leading to increasingly longer periodic cycles and finally deterministic chaos.

There are many models which demonstrate similar behavior (period-doubling bifurcations leading to deterministic chaos), for example, the Ricker model [33] 𝑥𝑛+1=𝑥𝑛𝑟exp1−𝑥𝑛(4.2) and the modification of the Hassel map [34]𝑓𝑟(𝑥)=𝑟𝑥1+(ğ‘Žğ‘¥)𝛾,𝑟>0,ğ‘Ž>0,𝛾>1.(4.3) Since in ecological systems chaotic dynamics is difficult to observe, it was suggested [35, 36] that perturbations can be a reason for this phenomenon. In fact, the perturbed Ricker model 𝑥𝑛+1=𝑥𝑛𝑟exp1−𝑥𝑛+𝜆,𝜆>0,(4.4) does not experience chaotic behavior for 𝑟 large enough; see [37] for some more details and relevant references. With the growth of 𝑟 the period-doubling route to chaos, which is characteristic to (4.4), will break down, giving rise to distinctive period-halving bifurcations and a stable two-cycle [35, 36] for 𝑟 large enough (and also for perturbation 𝜆<1, since for 𝜆>1 there is a stable equilibrium). Not all perturbed models demonstrate the same phenomenon for 𝑟 large. For example, the perturbed logistic model 𝑥𝑛+1=𝑟𝑥𝑛1−𝑥𝑛+𝜆,𝜆>0,(4.5) is chaotic for large 𝑟, whether 𝜆=0 or not.

Let us consider the perturbed version of model (2.1) 𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+ğ‘ğ‘’âˆ’ğ‘žğ‘¥(𝑛)+𝜆,𝑛≥0.(4.6) There are two parameters (𝑝 and ğ‘ž) in the model (4.6), and we will study bifurcation in both of them.

If we consider bifurcations of (2.1) in ğ‘ž, then there is again transition to deterministic chaos through the series of period-doubling bifurcations (see Figure 5(a)). The perturbed model (4.6) demonstrates a similar picture, up to a certain value of ğ‘ž; further, we observe period-doubling reversals leading to a stable equilibrium. This type of behavior can be rigorously justified.

Theorem 4.1. For fixed 𝑝>0, 𝛼∈(0,1), and 𝜆>0, there exists ğ‘ž0 such that for any ğ‘ž>ğ‘ž0 the positive equilibrium of (4.6) is globally asymptotically stable; moreover, the convergence is monotone.

Proof. We recall that the right-hand side of (4.6) 𝑔(𝑥)=(1−𝛼)𝑥(𝑛)+ğ‘ğ‘’âˆ’ğ‘žğ‘¥(𝑛)+𝜆 has a minimum point at 𝑥∗, where 𝑥∗=1ğ‘žî‚€lnğ‘ğ‘žî‚1−𝛼,limğ‘žâ†’âˆž1ğ‘žî‚€lnğ‘ğ‘žî‚1−𝛼=0.(4.7) Thus for fixed 𝑝, 𝛼, and 𝜆, there exists ğ‘ž0>0 such that 𝑥∗=1ğ‘žî‚€lnğ‘ğ‘žî‚1−𝛼<𝜆forğ‘ž>ğ‘ž0;(4.8)hence we also have 𝑓(𝑥∗)>𝜆>𝑥∗. This means that 𝑔(𝑥) is monotone increasing on the interval (𝑥∗,∞) including the equilibrium. By Remark 3.3 the only positive equilibrium of (4.6) is globally asymptotically stable; moreover, there is a monotone convergence to this equilibrium.

As far as bifurcations in ğ‘ž are considered, the perturbed equation (4.6) behaves like the perturbed version of the Ricker (4.2) and some other (4.3) population dynamics models: it demonstrates period-halving bifurcations and stable behavior for large ğ‘ž. If we follow bifurcations in 𝑝, then chaos is observed for any large 𝑝, similar to the perturbed logistic model (4.5).

5. Discussion and Open Problems

In the context of the model of red blood cells production, the mathematical properties of solutions can be interpreted in the following way.

Convergence to the equilibrium solution usually means healthy (normal) blood cell production, while sustainable oscillations correspond to a dynamic disease. In particular, if convergence to the equilibrium is monotone, then, in addition to normal blood cell production, parameter determination from data sampling is relatively easy compared to oscillatory cases. If the difference equation is considered as a numerical approximation of the continuous model, convergence to the equilibrium means that the numerical scheme is stable, in addition to the similar property of the continuous equation.

Deterministic chaos in the difference equation can mean one of the following two possibilities: first, blood cell production can be irregular (which occurs in practice, but was usually modeled by the addition of noise, or more general stochastic terms, or by introducing randomized parameters in the model), or, second, the numerical (discrete) approximation can create irregularity. In the case of the chaotic behavior, the range of parameters (which, as can be seen in Figures 4 and 5, can be separated from zero and bounded) are indicative for comparison with the normal (ideal) blood cells production. A small positive constant perturbation can correspond to a drug administration or other external influence stimulating blood cells production at a constant level.

The results of the present paper can be summarized as follows.(1)The model studied in the present paper can be considered as a discretization of the Lasota-Wazewska equation [3]. For the difference equation, we have presented solution bounds and strict nonoscillation conditions. Moreover, we have demonstrated that local asymptotic stability of the unique positive equilibrium is equivalent to its global asymptotic stability.(2)Bifurcations were investigated in both parameters 𝑝 and ğ‘ž involved in the equation. We also considered the linearly perturbed model. As 𝑝 grows, the bifurcation diagram of the perturbed equation does not differ qualitatively from the nonperturbed model (while demonstrating later bifurcations and transition to chaos than the original model). However, for large ğ‘ž the perturbed equation exhibits stable behavior. To the best of our knowledge, the chaotic behavior of (2.1) and (4.6) has never been studied before.

Finally, let us formulate some relevant open problems.(A)Consider bifurcations and chaos in a one-parameter version of (4.6) when neither 𝑝 nor ğ‘ž is fixed but there is a connecting relation 𝑝=𝑓(ğ‘ž). For example, similar to Theorem 4.1, the global stability can be easily justified for 𝑝=ğ‘žğ‘˜, 𝑘>0 large enough and 𝜆>0 in (4.6). Can a general theory be developed for some models similar to (2.1) which have a unique positive equilibrium and a unique critical point (which is a minimum rather than a maximum) and their linearly perturbed versions? Again, Theorem 4.1 demonstrates that, if the unique critical point tends to zero as the bifurcation parameter grows, then the perturbed model is stable for this parameter when it is large enough. Here we note that perturbed models with 𝜆<0 were investigated in [38] and encourage the study of models with negative constant perturbations which can describe immigration.(B)There are results for the analogue of (2.1) with periodic coefficients [27]. However, sharp oscillation and stability conditions are still an open problem.(C)Strictly speaking, the Lasota-Wazewska equation [3] involves a delay, so its discretization should be 𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+ğ‘ğ‘’âˆ’ğ‘žğ‘¥(𝑛−𝑘),𝑛≥0.(5.1) Prove that the local stability of the positive equilibrium of (5.1) implies its global stability.(D)As was mentioned in the beginning, the difference equation (2.1) was considered as a numerical discretization of the Lasota-Wazewska equation. A similar model for the Mackey-Glass equation was studied in [39], where the right-hand side is a unimodal function. If we consider the model corresponding to the Nicholson blowflies equation 𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+𝑝𝑥(𝑛−𝑘)ğ‘’âˆ’ğ‘žğ‘¥(𝑛−𝑘)(5.2) and its nondelay version 𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+𝑝𝑥(𝑛)ğ‘’âˆ’ğ‘žğ‘¥(𝑛),(5.3) then the right-hand side of (5.3), generally, has more than one critical point. However, it still has the only positive equilibrium, and it was demonstrated in [40] that local asymptotic stability implies global asymptotic stability of (5.3). Establish necessary and sufficient oscillation and global stability conditions for (5.2). Explore bifurcations and chaos for the perturbed version of (5.3).(E)The difference equation 𝐾𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+𝑝𝑥(𝑛)ln𝑥(𝑛)(5.4) can be considered as a discretization of the Gompertz differential equation. Obtain conditions on 𝛼,𝑝, and 𝐾 which guarantee that all solutions with 0<𝑥(0)<𝑀 for some 𝑀>0 remain positive. Deduce local stability conditions for the unique positive equilibrium and demonstrate that they imply global stability. Explore the delay version of this model 𝐾𝑥(𝑛+1)=(1−𝛼)𝑥(𝑛)+𝑝𝑥(𝑛)ln𝑥(𝑛−𝑘)(5.5) and deduce global stability results.


The first author was partially supported by the NSERC grant. The second author was partially supported by the Deanship of Scientific Research and the Research Centre in College of Science, King Saud University. The authors are grateful to anonymous referees for valuable comments and remarks which greatly improved the presentation of the paper.