A note on the Lasota discrete model for blood cell production

DISCRETE AND CONTINUOUS doi:10.3934/dcdsb.2019262 DYNAMICAL SYSTEMS SERIES B Volume 25, Number 2, February 2020 pp. 701–713 A NOTE ON THE LASOTA DISCRETE MODEL FOR BLOOD CELL PRODUCTION Eduardo Liz Departamento de Matem´atica Aplicada II, Universidade de Vigo 36310 Vigo, Spain Cristina Lois-Prados Instituto de Matem´aticas, Universidade de Santiago de Compostela Campus Vida, 15782 Santiago de Compostela, Spain Dedicated to Prof. Juan J. Nieto on the occasion of his 60th birthday Abstract. In an attempt to explain experimental evidence of chaotic oscillations in blood cell population, A. Lasota suggested in 1977 a discrete-time one-dimensional model for the production of blood cells, and he showed that this equation allows to model the behavior of blood cell population in many clinical cases. Our main aim in this note is to carry out a detailed study of Lasota’s equation, in particular revisiting the results in the original paper and showing new interesting phenomena. The considered equation is also suitable to model the dynamics of populations with discrete reproductive seasons, adult survivorship, overcompensating density dependence, and Allee effects. In this context, our results show the rich dynamics of this type of models and point out the subtle interplay between adult survivorship rates and strength of density dependence (including Allee effects). 1. Introduction. Based on experimental results, the Polish mathematician Andrzej Lasota proposed in 1977 [5] a discrete model for the production of red blood cells (erythrocytes). If xndenotes the number of cells at time n, and σxnis the number of cells destroyed in the time interval [n, n + 1], then the dynamics is governed by equation xn+1 −xn=−σxn+pn,(1.1) where pnis the number of cells which are produced in the bone marrow during the same period. The bone marrow has the ability to change production when the number of blood cells changes, so pn=p(xn). Based on experimental results, the form p(x)=(cx)γe−xwas proposed for the production function, thus leading to the difference equation xn+1 = (1 −σ)xn+ (cxn)γe−xn,(1.2) where σ∈(0,1) and c, γ are positive constants. As indicated by Lasota, equation (1.2) provides a flexible model that explains the behavior of the blood cell population in many clinical cases. For example, if γ > 1, then there is a critical value of the population size such that populations 2010 Mathematics Subject Classification. 39A10, 39A30, 92C37, 92D25. Key words and phrases. Discrete-time model, gamma model, Ricker map, stability, bifurcations, boundary collision, essential extinction. 701 702 EDUARDO LIZ AND CRISTINA LOIS-PRADOS below it cannot survive in the long term. In population dynamics, such effect is known as a strong Allee effect [4]. In the model under consideration, it means that if the number of cells is too small (for example, after a heavy hemorrhage), then the organism eventually dies. Cell populations above the critical value usually persist and they can approach a stable positive equilibrium (which is considered the normal behavior), but the number can oscillate when there is a disease affecting the blood cell production, and even behave in an erratic (aperiodic) way in the case of a severe disease. As the examples in [5] show, an increasing value of σleads to severe disease (because a large proportion of cells are destroyed). Even though this model assumes that the blood production is not a continuous process (a more accurate model seems to be a delay differential equation [5,16]), it is quite intuitive and displays a very rich dynamics. On the other hand, equation (1.2) can also be considered as a discrete model for iteroparous populations (that is, part of the population survives the reproductive season) whose recruitment function is defined by a gamma-Ricker map [8]. In this context, equation (1.2) is a flexible model for populations with discrete reproductive seasons, adult survivorship, overcompensating density dependence, and Allee effects. In this note, we show some interesting aspects of the dynamics of (1.2) that can help to understand the potential effects of increasing either the destruction rate of cells (in the erythropoietic model), or the adult mortality rate (in the population model). 2. Preliminary results. We first introduce the maps f: [0,∞)→[0,∞) and F: [0,∞)→[0,∞) defined by f(x)=(cx)γe−x, F(x) = (1 −σ)x+ (cx)γe−x= (1 −σ)x+f(x),(2.1) where σ∈(0,1), c, γ > 0. Notice that equation (1.2) can be written in the form xn+1 =F(xn). Some properties of the map fhave been recently proved in [8], and they will be useful in our analysis. Now, we state the intervals of monotonicity of Fdepending on the parameters. We prove that Fcan be either increasing or bimodal (see Figure 2.1). Theorem 2.1. The map Fdefined in (2.1) has the following properties: (a): F(0) = 0 and limx→∞ F(x) = ∞. (b): Fis differentiable in (0,∞),F0(0+)=1−σif γ > 1,F0(0+)=1−σ+cif γ= 1, and F0(0+) = ∞if 0< γ < 1. In any case, limx→∞ F0(x)=1−σ > 0. (c): Define H(γ) = √γ(γ+√γ)γ−1e−γ−√γ. (i) If cγH(γ)<1and 0< σ ≤σ∗:= 1 −cγH(γ), then Fis increasing in (0,∞). (ii) If cγH(γ)≥1, or cγH(γ)<1and σ > σ∗, then Fhas two positive critical points c1< c2. Moreover, Fattains a local maximum at c1and a local minimum at c2. Proof. The proof of items (a) and (b) is elementary, so we only provide the proof of (c). The first and second derivatives of Fin (0,∞) are F0(x) = (1 −σ) + cγxγ−1e−x(γ−x), F00(x) = cγe−xxγ−2(x2−2γx +γ2−γ). A NOTE ON THE LASOTA DISCRETE MODEL 703 0 1 2 3 4 5 6 0.0 0.5 1.0 1.5 2.0 σ < σ∗ σ > σ∗ Figure 2.1. Graphs of the map Ffor γ= 0.3, c= 0.6, and different values of σ:σ= 0.6< σ∗≈0.774 in red (Fis increasing), and σ= 0.9> σ∗in blue (Fhas two critical points). The dashed line is the graph of y=x. We distinguish two cases according to the zeros of F00 on (0,∞) (see Figure 2.2). (a) x (b) x 0 1 2 3 4 5 6 7 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 7 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Figure 2.2. Graphs of the map F0(solid, blue) and the line y= 0 (dashed, black), with c=1. (a): γ= 0.25 <1, σ= 0.8> σ∗≈ 0.707; (b): γ= 2 >1, σ= 0.9> σ∗≈0.841. In both cases, there are two intersection points, that determine the critical points of F. If γ≤1, then F0has a unique critical point at d+=γ+√γ; moreover, F0is decreasing on (0, d+) and increasing on (d+,∞). If γ > 1, then F00(d+) = F00(d−) = 0, where d+=γ+√γ > d−=γ−√γ > 0. In this case, F0(0+)>0, F0is increasing on (0, d−)∪(d+,∞), and decreasing on (d−, d+). In both cases, it is clear that F0(x)≥0 for all x > 0 if and only if F0(d+)≥0, and F0has two positive zeros c1< c2if F0(d+)<0. Moreover, it follows from the sign of F0that c1corresponds to a local maximum of F, and c2to a local minimum. Finally, it is easy to check that condition F0(γ+√γ)≥0 is equivalent to σ≤ 1−cγH(γ). 3. Fixed points and stability. We first address the case 0 < γ ≤1, for which the map Fcan have at most one positive equilibrium. 704 EDUARDO LIZ AND CRISTINA LOIS-PRADOS Theorem 3.1. Assume that 0< γ ≤1. (a): If γ < 1, then Fhas a unique positive fixed point p. Moreover, F(x)> x if 0< x < p, and F(x)< x if x>p. (b): If γ= 1, then Fhas a unique positive fixed point pif c>σand no positive fixed points if c≤σ. In the former case, F(x)> x if 0< x < p, and F(x)< x if x>p. (c): Assume that 0< γ < 1, or γ= 1 and c>σ. Then the unique positive equilibrium pis locally asymptotically stable for the Lasota equation (1.2) if the following condition holds: cγ< G(γ, σ):=σγ+2 σ−11−γ eγ+2 σ−1.(3.1) If cγ> G(γ, σ), then pis unstable. Proof. We only prove the case γ < 1, since the proof for γ= 1 is very similar. The fixed point p > 0 satisfies the identity f1(p):=σp1−γ=cγe−p=:f2(p). The existence and uniqueness of pfollows from the fact that f1(x) = σx1−γis increasing, with f1(0) = 0, and f2(x) = cγe−xis decreasing, with f2(0) >0, f2(∞) = 0. A simple graphical analysis shows that F(x)> x if 0 < x < p, and F(x)< x if x>p. The derivative of Fat the fixed point pis F0(p)=1−σ+γcγpγ−1e−p−cγpγe−p= 1 + σ(γ−1−p).(3.2) Since γ < 1, it is clear that F0(p)<1. On the other hand, it follows from (3.2) that F0(p)>−1⇐⇒ p<γ+2 σ−1. Using the statement in (a), it is clear that the last inequality holds if and only if Fγ+2 σ−1< γ +2 σ−1, and this inequality is equivalent to (3.1). Remark 3.2. We notice that equation (1.2) is uniformly permanent if either γ < 1 or γ= 1 and σ < c. This means that there exists a compact interval [A, B] such that 0 < A ≤lim infn→∞ Fn(x)≤lim supn→∞ Fn(x)≤B, for all x > 0. Indeed, if Fis increasing then the unique positive equilibrium pis a global attractor, and we can choose A=B=p; if Fis bimodal, then elementary arguments show that we can choose A=F(c2), B=F(c1), where c1,c2are the critical points of F (0 < c1< c2). Next, if γ > 1, then Fcan have 0, 1, or 2 positive fixed points. This case is studied in our next result. Theorem 3.3. Assume that γ > 1, and define c∗:= σ(γ−1)1−γeγ−11/γ .(3.3) Then the map Fdefined in (2.1) has two positive fixed points if c > c∗, one positive fixed point if c=c∗, and no positive fixed points if c < c∗. If Fhas only one fixed point q, then q=γ−1,F0(q) = 1 and qis semi-stable. If Fhas two positive fixed A NOTE ON THE LASOTA DISCRETE MODEL 705 points q, p, then q < γ −1< p,F0(q)>1, and F0(p)<1. The fixed point qis unstable, while pis asymptotically stable if (3.1) holds. If c<c∗, then all solutions of (1.2) converge to 0. Proof. The positive fixed points of Fare the solutions of equation f1(x) = f2(x), where f1(x) = σx1−γand f2(x) = cγe−x. Since f1(0+) = ∞,f2(0) = cγ,f1 and f2are decreasing and convex, and limx→∞(f1(x)/f2(x)) = ∞, it follows from elementary arguments that Fcan have at most two positive fixed points. Assume that there is at least one positive fixed point of F, and denote by q the smallest one. It follows from the properties of f1and f2that f0 1(q)≤f0 2(q). This inequality, together with the equality f1(q) = f2(q), leads to q≤γ−1. If q=γ−1, then f0 1(q) = f0 2(q) and γ−1 is the unique positive fixed point of F. However, if q < γ −1, then f0 1(q)< f0 2(q). Since limx→∞(f1(x)/f2(x)) = ∞, there exists another fixed point p>q. Moreover, f0 1(p)> f0 2(p), which is equivalent to p>γ−1. Thus: •there is a unique positive fixed point of Fif and only if q=γ−1, that is, if f1(γ−1) = f2(γ−1). This equality is equivalent to c=c∗; •there are two positive fixed points q, p of Fif and only if f1(γ−1) < f2(γ−1), which is equivalent to c > c∗; •Fdoes not have positive fixed points if c < c∗. In the first case, it is easy to check that solutions of (1.2) starting at an initial condition x0∈(0, q) converge to zero. If F(x)≥qfor all x > q, then all solutions starting at [q, ∞) converge to q; if there are points x>qsuch that F(x)< q, then the immediate basin of attraction of qis [q, r], where ris the smallest point in F−1(q)\{q}. Therefore, qis semistable. In the second case, we have that F0(p)<1< F0(q). Thus, qis unstable. If (3.1) holds, then F0(p)>−1, implying that the equilibrium pis asymptotically stable. Finally, if c < c∗, then F(x)< x for all x > 0 and therefore all solutions of (1.2) converge to 0. Figure 3.1 shows the three different possibilities considered in the statement of Theorem 3.3. It is well known that a discrete dynamical system generated by a map with two critical points can lead to the coexistence of several attractors. However, if γ < 1 or γ= 1 and σ < c, we will prove that the unique positive fixed point pis a global attractor of all positive solutions of (1.2) when it is asymptotically stable. Moreover, we get that the parameter region of global asymptotic stability includes the nonhyperbolic case, thus extending condition (3.1) to cγ≤G(γ, σ). To prove this result, we need the following generalization of Theorem 1 in [10]. Theorem 3.4. Assume that 0< σ < 1and f: [0,∞)→[0,∞)satisfies the following conditions: (A1): g(x) = (1/σ)f(x)has a unique positive fixed point p > 0,g(0) = 0, and g0(0+)>1(g0(0+)can be ∞). (A2): fhas a unique critical point z; moreover, f0(x)>0for all x∈(0, z)and f0(x)<0for all x>z. (A3): (Sf)(x)<0for all x>z, where (Sf)(x) = f000(x) f0(x)−3 2f00(x) f0(x)2 706 EDUARDO LIZ AND CRISTINA LOIS-PRADOS is the Schwarzian derivative of f. (A4): f00(x)<0for all x∈(0, z). Then, the unique positive equilibrium pof equation xn+1 = (1 −σ)xn+f(xn) (3.4) is globally asymptotically stable if 1−σ+f0(p)≥ −1.(3.5) If (3.5) does not hold, then pis unstable. 0 5 10 15 20 25 0 5 10 15 20 25 c=c∗ c < c∗ c > c∗ Figure 3.1. Graphs of the map Ffor γ= 8, σ= 0.8, and different values of c:c=c∗≈0.425 in blue (one positive fixed point), c= 0.4< c∗in red (no positive fixed points), and c= 0.5> c∗in black (two positive fixed points). The dashed line is the graph of y=x. Proof. Equation (3.4) can be written in the form of equation (4) in [10]: xn+1 =αxn+ (1 −α)g(xn),(3.6) with α= 1 −σand g(x) = (1/σ)f(x). It is clear that conditions (A2)–(A4) hold for gbecause g0(x) = (1/σ)f0(x), g00(x) = (1/σ)f00(x), and (Sg)(x)=(Sf)(x). There are two relevant differences with Theorem 1 in [10] . On the one hand, condition (A3) there required (Sf)(x)<0 for all x6=z. However, a simple inspection of the proof shows that the less restrictive condition (Sf)(x)<0 for all x>zis enough to get the result. On the other hand, the restriction z < p is required in [10, Theorem 1]. In case z≥p, the map φ(x) = (1 −σ)x+f(x) defining the right-hand side of (3.4) satisfies φ0(x)>0 for 0 < x ≤p,φ(x)> x for x < p, and 0 < φ(x)< x for x>p. Hence, simple arguments show that pis globally asymptotically stable (see, e.g., [3, Lemma 1]). Now we are in a position to prove the global stability result for (1.2). Corollary 3.5. Assume that 0< γ < 1, or γ= 1 and c > σ. Then, the unique positive equilibrium pfor the Lasota equation (1.2) is globally asymptotically stable if and only if cγ≤G(γ, σ), where G(γ, σ)is defined in (3.1). Proof. Consider the map f(x)=(cx)γe−xdefined in (2.1). It is clear from the proof of Theorem 3.1 that condition cγ≤G(γ, σ) is equivalent to (3.5). Thus, in order to use Theorem 3.4, we only need to check conditions (A1)–(A4). A NOTE ON THE LASOTA DISCRETE MODEL 707 Conditions (A1) and (A2) follow from elementary arguments, and (A3) was proved in [8, Proposition 2]. Hence, it remains to prove (A4). From the proof of Theorem 2.1, we know that f00(x)<0 for all x∈(0, γ +√γ). Since the unique critical point of fis z=γ < γ +√γ,(A4) follows. 4. Dynamics and bifurcations. In this section, we combine the rigorous analysis from the previous sections with numerical bifurcation diagrams to show the potential rich dynamics of the solutions to the Lasota equation (1.2) and the role of the parameters γand σ. For parameter c, we fix the value c= 0.47 already used by Lasota in [5]. Although his study is restricted to γ= 8 and three different values of σ(σ= 0.1,0.4,0.8), we will consider more cases to give an idea of the more relevant dynamical phenomena. In particular, we will revisit the cases studied by Lasota. 0 2 4 6 8 10 12 0.0 0.2 0.4 0.6 0.8 1.0 σ γ c Extinction Bistability Oscillations and Essential Extinction Global Stability Figure 4.1. Main bifurcation boundaries and regions with different dynamical behavior for equation (1.2) with c= 0.47, in the parameter plane (γ, σ). The two solid lines represent the extinction boundary (red color) and the stability boundary of the largest positive equilibrium (blue color). The vertical dashed line γ= 1 (from σ= 0 to σ=c= 0.47) is the border between global stability of the unique positive equilibrium and a bistability region, in which both the largest positive equilibrium and the extinction equilibrium are asymptotically stable. We begin by displaying in Figure 4.1 a bifurcation diagram for c= 0.47, which shows in the parameter plane (γ, σ) the stability and extinction boundaries provided in Section 3. In this case, condition cγ≤G(γ, σ) holds for all γ∈(0,1], and therefore Corollary 3.5 implies that the unique positive equilibrium is globally asymptotically stable in this parameter range (if γ= 1, this equilibrium exists only for σ < 0.47 = c). The red solid line is given by equation (3.3) when γ > 1. According to Theorem 3.3, this line corresponds to a tangent bifurcation and represents the boundary between extinction and bistability. In the bistable case, there are solutions that converge to the extinction equilibrium and others that converge to the largest positive 708 EDUARDO LIZ AND CRISTINA LOIS-PRADOS equilibrium, so permanence depends on the initial condition. For σ > c, the red line γ= 1 represents the boundary between global stability of the positive equilibrium and extinction. The blue solid line represents the boundary of local asymptotic stability of p, so the region to the right of this line corresponds to more complicated dynamics, which includes periodic attractors, chaotic attractors, and essential extinction, as we show below in some numerical bifurcation diagrams. We recall that essential extinction means that almost all solutions (in a sense of Lebesgue measure) converge to zero [14,15]. Essential extinction usually occurs after a boundary collision between the basins of attraction of a chaotic attractor and the extinction equilibrium 0; however, we report below different transitions from bistability to essential extinction. In the remainder of this section, we show some relevant numerical bifurcation diagrams to illustrate the rich behavior of equation (1.2). The value of cis always set to 0.47. In some cases (Figures 4.2,4.3), we fix a value of γand use σas a bifurcation parameter, while in Figure 4.4 we do the opposite. 4.1. Bubbles. A first interesting phenomenon is the existence of the so-called bubbles, which are characteristic in bifurcation diagrams of population models where some adults survive more than one reproduction period. A primary bubble in the sense of [11, Definition 3] is shown in Figure 4.2: the largest positive equilibrium loses its asymptotic stability in a period-doubling bifurcation as the destruction rate σis increased, but stability is regained at a larger value of σ, after a periodhalving bifurcation occurs. This phenomenon shows that increasing the mortality rate can be either stabilizing or destabilizing. More complex bubbling effects can be observed in Figure 4.3 (c). In this case, the equation can reach a chaotic regime inside the bubble (see [11] for more details and related references). 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 population size, xn destruction rate, σ Figure 4.2. Bifurcation diagram showing a bubble for equation (1.2) with c= 0.47, γ= 7.65, using σas the bifurcation parameter. Black dashed lines correspond to unstable equilibria. 4.2. Hydra effect. It is easy to check that the largest equilibrium of (1.2) decreases as σis increased. However, the average population size of the nontrivial attractor can increase with σwhen the largest equilibrium loses its asymptotic stability. This phenomenon can be observed in Figure 4.3 (b), where the average of the A NOTE ON THE LASOTA DISCRETE MODEL 709 2-periodic attractor is represented by a solid red line. The phenomenon of a population increasing in response to an increase in its per-capita mortality rate is known as the hydra effect [1]. A formal notion of hydra effect suitable for the context of this paper is given in [11, Definition 2]. population size, xn population size, xn destruction rate, σ (a) destruction rate, σ (b) 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 population size, xn population size, xn destruction rate, σ (c) destruction rate, σ (d) 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 25 30  0 5 10 15 20 25 30 0.825 0.83 0.835 Figure 4.3. Bifurcation diagrams for equation (1.2) with c= 0.47, using σas the bifurcation parameter. Black dashed lines correspond to unstable equilibria. For more details, see the text. (a): γ= 7; (b): γ= 8; (c): γ= 8.5 and σ∈(0,1); (d): magnification for γ= 8.5. 4.3. Sudden collapses. A typical feature of population models with Allee effects, where population cannot survive in the long term if its abundance is below a critical size, is the possibility of sudden collapses. Roughly speaking, this means that there are critical values of some relevant parameter such that extinction occurs if the parameter goes beyond this critical value, while populations can persist at densities bounded away from zero for values just below the critical one. Sudden collapses typically occur when the basins of attraction of the zero equilibrium and another attractor collide (boundary collisions). See, e.g., [6,14,15]. There are two main mechanisms leading to sudden collapses. The simplest one is a tangent bifurcation, which leads from bistability to extinction when the two positive equilibria of (1.2) collide and then disappear (see Figure 3.1). A bifurcation