\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 140, pp. 1--20.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/140\hfil Age-structured model of hematopoiesis dynamics] {Age-structured model of hematopoiesis dynamics with growth factor-dependent coefficients} \author[M. Adimy, Y. Bourfia, M. L. Hbid, C. Marquet \hfil EJDE-2016/140\hfilneg] {Mostafa Adimy, Youssef Bourfia, \\ My Lhassan Hbid, Catherine Marquet} \address{Mostafa Adimy \newline INRIA, Universit\'e Lyon, Universit\'e Lyon 1, Institut Camille Jordan, 43 Bd. du 11 novembre 1918, F-69200 Villeurbanne Cedex, France} \email{mostafa.adimy@inria.fr} \address{Youssef Bourfia \newline Universit\'e Cadi Ayyad, Facult\'e des Sciences Semlalia, D\'epartement de Math\'ematiques, UMMISCO-Marrakech (IRD-UPMC), BP. 2390, Marrakesh, Morocco. \newline Universit\'e Pierre et Marie CURIE, Laboratoire Jacques-Louis Lions, CNRS UMR 7598, place Jussieu 75005 Paris, France} \email{bourfiayoussef@gmail.com} \address{My Lhassan Hbid \newline Universit\'e Cadi Ayyad, Facult\'e des Sciences Semlalia, D\'{e}partement de Math\'ematiques, UMMISCO-Marrakech (IRD-UPMC), BP. 2390, Marrakesh, Morocco} \email{hbid@uca.ma} \address{Catherine Marquet \newline Universit\'e de Pau, Laboratoire de Math\'ematiques Appliqu\'ees, CNRS UMR 5142, Avenue de l'universit\'e, 64000 Pau, France} \email{catherine.marquet@univ-pau.fr} \thanks{Submitted February 9, 2016. Published June 10, 2016.} \subjclass[2010]{34D20, 34D23, 34K06, 92C37} \keywords{Age-structured partial differential equation; delay differential system; \hfill\break\indent Lyapunov function; characteristic equation; cell dynamic} \begin{abstract} We propose and analyze an age-structured partial differential model for hematopoietic stem cell dynamics, in which proliferation, differentiation and apoptosis are regulated by growth factor concentrations. By integrating the age-structured system over the age and using the characteristics method, we reduce it to a delay differential system. We investigate the existence and stability of the steady states of the reduced delay differential system. By constructing a Lyapunov function, the trivial steady state, describing cell's dying out, is proven to be globally asymptotically stable when it is the only equilibrium of the system. The asymptotic stability of the positive steady state, the most biologically meaningful one, is analyzed using the characteristic equation. This study may be helpful in understanding the uncontrolled proliferation of blood cells in some hematological disorders. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \allowdisplaybreaks \section{Introduction}\label{S1} Hematopoiesis is the physiological process that ensures the production and regulation of blood cells. It involves a pool of undifferentiated and self-renewing cells called hematopoietic stem cells (HSCs), located in the bone marrow, from which arises all differentiated blood cell lineages (red blood cells, white cells and platelets). Proliferation, differentiation and apoptosis are processes occurring during hema\-topoiesis and are all mediated by a wide range of hormone-like molecules called growth factors. The growth factors play an activator or inhibitor role in this process and they act on every cell compartment: primitive stem cells, progenitors and precursors. Their role to maintain homeostasis of blood cells is essential. The production of red blood cells (erythropoiesis) and platelets (megakaryopoiesis) seems to be regulated by specific growth factors whereas white blood cell production (leukopoiesis) is more complicated and less clearly understood. For the red blood cells, the erythropoietin (EPO) helps to regulate erythrocyte production (Adamson \cite{AdamsonAJM1996}). A decrease in mature red blood cell count leads to a decrease in tissue $\textit{p}0_2$ levels, which in turn increases the production of EPO by the kidneys and controls erythropoiesis. For the platelets, it seems that their production and regulation are controlled by feedback mechanisms involving specific cytokines such as thrombopoietin (TPO). However, it has been shown that the cytokine TPO affects other cell lines as well (Tanimukai et al \cite{TanimukaiKimSakAlEH1997}), which means that the three lines are probably not fully independent, and there is a feedback control from mature cells to HSCs. Regulation of the multiple fates of HSCs, including quiescence, self-renewal, differentiation and apoptosis, requires the cooperative actions of several growth factors that bind to receptors on these cells. Many of the important players in this regulation have been identified (Tanimukai et al \cite{TanimukaiKimSakAlEH1997}). Due to the number of divisions, and the quantity of cells and cytokines involved in hematopoiesis, issues may arise at different cellular levels and sometimes result in disorders affecting blood cells. Among a wide variety of disorders affecting blood cells, myeloproliferative diseases are of great interest. They are characterized by a group of conditions that cause blood cells to grow abnormally. They include chronic myelogenous leukemia, a cancer of white blood cells. In some cases, chronic myelogenous leukemia exhibits periodic oscillations in all blood cell counts (see \cite{PujoBernardMackeySJADS2005}). Myeloproliferative disorders usually originate from the HSC compartment: an uncontrolled proliferation in the HSC compartment can perturb the entire system and leads to a quick or slow proliferation. A low blood counts (white cell count, red cell count, or platelet count) can be associated with many diseases and conditions that cause the body to have too few blood cells. It can be associated to a bone marrow failure, consecutive to disease of another organ (for example, liver or kidney), or secondary to treatment with some drugs (for example, chemotherapy drugs). Mathematical modeling of hematopoiesis dynamics has been the focus of a large panel of researchers over the last four decades, with attempts to improve the understanding of the complex mechanisms regulating HSC functions, throughout the course of normal and pathological hematopoiesis. One of the earliest mathematical models that shed some light on this process was proposed by Mackey \cite{Mackey1978} in 1978 inspired by the work of Lajtha \cite{Lajtha1959}, and Burns and Tannock \cite{BurnsTannock1970}. Mackey's model is a system of two delay differential equations describing the evolution of the HSC population divided into proliferating and quiescent cells (also called resting cells). This model has been studied, analyzed and applied to hematological diseases by many authors (see for instance, \cite{AdimyCrausteRuanSIAP2005, MackeyRudnickiJMB1994, MackeyRudnickiJMB1999, PujoBernardMackeySJADS2005, PujoMackeyCRB2004, ReyMackeyCAMQ1993}). We refer the reader interested in this topic, in addition to the previous articles, to the review papers by Adimy and Crauste \cite{AdimyCrausteMMNP2012}, Haurie et al \cite{HaurieDaleMackeyBLOOD1998}, Mackey et al \cite{MackeyHaurieBelairBOOK2003}, and the references therein. Mathematical models describing the action of growth factors on the hematopoiesis process have been proposed by B\'elair et al in 1995 \cite{BelairMackeyMahaffyMB1995}, and Mahaffy et al in 1998 \cite{MahaffyBelairMackeyJTB1998}. They considered an age-structured model of HSC dynamics, coupled with a differential equation to describe the action of a growth factor on the reintroduction rate from the resting phase to the proliferating one. In 2006, Adimy et al \cite{AdimyCrausteRuanBMB2006} proposed a system of three delay differential equations describing the production of blood cells under the action of growth factors assumed to act on the rate of reintroduction into the proliferating phase. Adimy and Crauste considered and analyzed two models of hematopoiesis dynamics with: the influence of growth factors on HSC apoptosis \cite{AdimyCrausteDCDSB2007}, and the action of growth factors on the apoptosis rate as well as on the reintroduction rate into the proliferating phase \cite{AdimyCrausteMCM2009}. In this paper, we consider the influence of growth factors on the apoptosis rate, on the differentiation rates (of the proliferating and quiescent cells), as well as on the reintroduction rate into the proliferating phase (see Figure \ref{Fig1}). The resulting system is composed by three age-structured partial differential equations for the different compartments of cell population, coupled with a system of four differential equations to describe the action of growth factors on different parameters of the system. To our knowledge this model has never been considered in hematopoiesis dynamics. The paper is organized as follows. In section \ref{S2}, we provide some biological background leading to an age-structured partial differential model for HSC dynamics. In section \ref{S3}, we use the method of characteristics to reduce the model to a system of delay differential equations. In section \ref{S4}, we establish some proprieties of the solutions such as positivity and boundedness. In section \ref{S5}, we investigate the existence of steady states. In section \ref{S6}, we prove the global asymptotic stability of the trivial steady state using a Lyapunov function. In section \ref{S7}, we linearize the delay system about each steady state and we deduce the delay-dependent characteristic equation. Then, we obtain the local asymptotic stability of the positive steady state. \section{Age-structured partial differential model}\label{S2} We consider two cell populations, HSC population (in the bone marrow) and mature blood cell population (in the bloodstream), for instance red blood cells. The HSC population is divided into proliferating and quiescent cells. Proliferating cells are the ones performing the cell division (growth, DNA synthesis and mitosis). Quiescent (or resting) HSCs are actually in a quiescent phase ($G_0$-phase). HSCs generate cells that undergo terminal differentiation resulting in mature circulating blood cells. Mature blood cells control the HSC population through growth factors. We denote respectively by $n(t,a)$, $p(t,a)$ and $m(t,a)$ the cell population densities of quiescent HSCs, proliferating HSCs and mature cells, with age $a \geq 0$ at time $t \geq 0$. The age represents the time spent by a cell in one of the three compartments. A schematic representation of this model is given in Figure \ref{Fig1}. Details of the modeling are presented hereafter. Quiescent cells are assumed to die with a constant rate $0 \leq \delta \leq 1$, and they can be introduced into the proliferating phase with a rate $\beta$ in order to divide. We suppose that $\beta$ depends upon a growth factor concentration $E_1$, that stimulates the proliferative capacity of HSCs: the more growth factor, the more proliferation of HSCs. Hence the feedback induced by the growth factor $E_1$ is positive, and the function $\beta$ is supposed to be increasing, with $\beta(0)=0$. As soon as a cell enters the proliferating phase, it is committed to divide a time $\tau \geq 0$ later. We assume that the duration of the proliferating phase is the same for all cells, so $\tau$ is constant, and describes an average duration of the cell cycle. The population of proliferating cells is controlled by apoptosis $\gamma \geq 0$, which is a programmed cell death that eliminates deficient cells and also maintains the homeostatic state of cell population. We assume that the apoptosis rate $\gamma$ depends upon the concentration of growth factor $E_2$ (for example, EPO, see \cite{KouryBondurant1990}). Since an increase of the growth factor concentration $E_2$ leads to a decrease of the apoptosis rate, we assume that $\gamma$ is a decreasing function of $E_2$ and $\lim_{E_2 \to +\infty}\gamma(E_2)=0$. The portion of quiescent cells that differentiate to mature cells is denoted by $K_N \geq 0$ which, we assume, depends upon a growth factor concentration denoted $E_3$. Since an increase of growth factor $E_3$ leads to an increase of the differentiation, we suppose that $E_3 \mapsto K_N(E_3)$ is an increasing function. Here, we only consider one kind of mature cells, for instance red blood cells. Then, we can consider that $1-(\delta + \beta + K_N) \geq 0$, the remainder of quiescent cells, differentiate to other cell lineages (for instance, white blood cells and platelets). At the end of the proliferating phase, each cell divides into two daughter cells. The daughter cells can either differentiate and enter the mature phase or stay in HSC compartment and enter to the $G_0$-phase. We assume that the part of daughter cells $\alpha \geq 0$ that stay in HSC compartment is constant. This is important because HSCs could maintain their characteristic properties of self-renewal and lack of differentiation could provide an unlimited source of cells to maintain the homeostasis. The portion of daughter cells entering the mature phase is denoted by $K_P \geq 0$ which, we assume, depends upon a growth factor concentration denoted $E_4$. As for quiescent cells, we suppose that $E_4 \mapsto K_P(E_4)$ is an increasing function and that the portion $1-(\alpha + K_P) \geq 0$ of daughter cells differentiate to other cell lineages. We suppose that the mature cells die with a constant rate $\mu \geq 0$. All the growth factor concentrations $E_1$, $E_2$, $E_3$ and $E_4$ are controlled by the mature cells through functions $f_i$, $i=1,2,3,4$ acting as negative feedbacks of the mature blood cells on the production of growth factors (see Figure \ref{Fig1}). \begin{figure}[htb] \begin{center} \includegraphics[width=0.9\textwidth]{fig1} % scheme.pdf \end{center} \caption{Schematic representation of HSC dynamics. Solid arrows represent the mechanisms taken into account: differentiation, cell division, reintroduction into the proliferating phase, apoptosis and natural death. The dependency of the parameters upon growth factors are represented by dashed lines. The dash-dotted line represents the feedback control from mature cells to the growth factors.} \label{Fig1} \end{figure} The densities $n(t,a)$, $p(t,a)$ and $m(t,a)$ satisfy, for $t>0$, the system \begin{equation} \begin{gathered} {\frac{\partial n}{\partial t}} + {\frac{\partial n}{\partial a}} = -(\delta+\beta(E_1(t)))n(t,a), \quad a > 0, \\ {\frac{\partial p}{\partial t}} + {\frac{\partial p}{\partial a}} = -\gamma(E_2(t))p(t,a), \quad 0 < a < \tau, \\ {\frac{\partial m}{\partial t}} + {\frac{\partial m}{\partial a}} = -\mu m(t,a), \quad a > 0. \end{gathered} \label{eq.mdl.1} \end{equation} System \eqref{eq.mdl.1} is completed by boundary conditions (for $a = 0$), that describe the flux of cells entering each phase, and by initial conditions (for $t=0$). Then the boundary conditions of \eqref{eq.mdl.1} are, for $t>0$, \begin{equation} \begin{gathered} n(t,0) = 2\alpha p(t,\tau),\\ p(t,0) = \beta(E_1(t))N(t),\\ m(t,0) = K_N(E_3(t))N(t)+2K_P(E_4(t))p(t,\tau), \end{gathered}\label{cond.l.1} \end{equation} where \begin{equation*} N(t)=\int_0^{+\infty}{n(t,a)da}, \quad P(t)=\int_0^{\tau} {p(t,a)da}, \quad M(t)=\int_0^{+\infty} {m(t,a)da}, \end{equation*} and $E_i(t)$, $i=1, 2, 3, 4$, are growth factor concentrations. Initial conditions of \eqref{eq.mdl.1} are given by nonnegative $L^1$-functions $n_0$, $p_0$ and $m_0$, such that \begin{equation} \begin{gathered} n(0,a) = n_0(a), \quad m(0,a) =m_0(a), \quad \text{for } a \geq 0, \\ p(0,a)=p_0(a), \quad \text{for } a \in [0,\tau]. \end{gathered} \label{cond.i} \end{equation} In addition, we assume that \begin{equation*} \lim_{a \to+ \infty } n( {t,a} ) = \lim_{a \to +\infty } m( {t,a} ) = 0, \quad \text{for } t \geq 0. \end{equation*} The growth factor concentrations $E_i(t)$, follow the evolution equations \begin{equation}\label{eq.factors} E'_i(t) = -k_iE_i(t) + f_i(M(t)), \end{equation} where the coefficients $k_i > 0$ are the degradation rates of the growth factors $E_i$. We assume that the functions $M \mapsto f_i(M)$ are positive, decreasing and satisfy $\lim_{M \to +\infty } f_i(M) = 0$. \section{Reduction to a delay differential system}\label{S3} The age-structured model \eqref{eq.mdl.1}-\eqref{cond.l.1}-\eqref{cond.i}-\eqref{eq.factors} can be reduced to a delay differential system. The method of characteristics implies, for $t > 0$ and $a \in (0,\tau)$, that \begin{equation}\label{solutions.de.P} p(t,a)=\begin{cases} p_0(a-t)\exp\Big(-\int_0^{t}{\gamma(E_2(s))ds}\Big), & \text{if } 0\tau. \end{cases} \end{aligned} \end{equation} Integrating the last equation of \eqref{eq.mdl.1} with respect to the age variable, we obtain \begin{equation}\label{intg3} M'(t) = - \mu M( {t} ) +m( {t,0} ). \end{equation} Then, using the last boundary condition of \eqref{cond.l.1}, we obtain \begin{equation*} M'(t) = - \mu M( {t} ) + K_N(E_3(t))N(t) + 2K_P(E_4(t))p(t,\tau). \end{equation*} We conclude that \begin{equation}\label{eq.imp.suit3} \begin{aligned} M'(t) & = - \mu M( {t} ) + K_N(E_3(t))N(t) + 2K_P(E_4(t)) \\ &\quad\times \begin{cases} p_0(a-t)\exp(-\int_0^{t}{\gamma(E_2(s))ds}), & \text{if } t<\tau, \\ \beta(E_1(t-\tau)) N(t-\tau) \exp ( - \int_{t-\tau}^{t} { \gamma(E_2(s))\,ds}), & \text{if }t>\tau. \end{cases} \end{aligned} \end{equation} Note that, for $t>\tau$, we have \begin{equation*} P(t)=\int_0^{\tau}\beta(E_1(t-a))N(t-a)\exp\Big(-\int_{t-a}^{t}{\gamma(E_2(s))ds}\Big)da. \end{equation*} Then, the asymptotic behavior of $P$ is related to $E_1$, $E_2$ and $N$. On the other hand, $N$, $M$, and $E_i$, do not depend on $P$, then, we can focus on the study of the solutions $(N,M,E_i)$. One can notice that, on the interval $[0,\tau]$ the functions $(N,M,E_i)$ satisfy a non-autonomous ordinary differential system, and for $t>\tau$, they satisfy the delay differential system \begin{equation}\label{sys.imp.suit} \begin{aligned} N'(t) &= - (\delta + \beta(E_1(t))) N( {t} ) \\ &\quad + 2\alpha\beta(E_1(t-\tau)) N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E_2(s))\,ds}\Big),\\ M'(t) &= - \mu M( {t} ) + K_N(E_3(t))N(t) \\ &\quad + 2K_P(E_4(t))\beta(E_1(t-\tau)) N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E_2(s))\,ds}\Big),\\ E'_i(t) &= -k_iE_i(t) + f_i(M(t)), \end{aligned} \end{equation} with initial conditions solutions of the ordinary differential system \eqref{eq.factors}-\eqref{eq.imp.suit}-\eqref{eq.imp.suit3} defined on the interval $[0,\tau]$. For each continuous initial condition, the system \eqref{sys.imp.suit} has a unique solution, defined for $t>\tau$ (see Hale and Verduyn Lunel \cite{Hale1993S}). From now on, we make a translation of the initial conditions so as to define them on the interval $[-\tau,0]$, as it can be found in Hale and Verduyn Lunel \cite{Hale1993S}. \section{Positivity and boundedness of solutions}\label{S4} We focus on the positivity and boundedness properties of the solutions $(N,M,E_i)$ of system \eqref{sys.imp.suit}. The following result states that all solutions of system \eqref{sys.imp.suit} are nonnegative, provided that initial conditions are nonnegative. \begin{proposition} The solutions of system \eqref{sys.imp.suit} associated with nonnegative initial conditions are nonnegative. \end{proposition} \begin{proof} Let $(N(t),M(t),E_i(t))$ be a solution of \eqref{sys.imp.suit}. Firstly, we check that $N$ is nonnegative. Assume that there exist $t_0>0$ and $\epsilon \in (0,\tau)$ such that $N(t)>0$, for $0 < t < t_0$, $N(t_0) = 0$ and $N(t)<0$, for $t\in (t_0, t_0+\epsilon)$. Let $t\in (t_0, t_0+\epsilon)$. It follows from \eqref{sys.imp.suit}, that \begin{equation*} N'(t_0) = 2\alpha\beta(E_1(t_0-\tau)) N(t_0-\tau) \exp \Big( - \int_{t_0-\tau}^{t_0} { \gamma(E_2(s))\,ds}\Big) >0. \end{equation*} This gives a contradiction. Consequently, $N(t)$ is nonnegative for $t\geq0$. Using a similar reasoning, we prove that $M(t)$ is nonnegative. Finally, the positivity of $E_i(t)$ follows from the fact that $f_i$ is positive. \end{proof} We now concentrate on the boundedness properties of the solutions of system \eqref{sys.imp.suit}. We start by proving the following lemma. \begin{lemma}\label{L.4.2} The solution $E_i(t)$ of \eqref{sys.imp.suit} is strictly decreasing as long as $E_i(t)>f_i(0)/k_i$, and either: \begin{itemize} \item[(i)] $E_i(t) > f_i(0)/k_i$ for all $t \geq 0$ and then, $\lim_{t \to +\infty}E_i(t)=f_i(0)/k_i$, or \item[(ii)] there exists $\overline{t}_i \geq 0$ such that $E_i(\overline{t}_i) = f_i(0)/k_i$ and then, $E_i(t) \leq f_i(0)/k_i$, for all $t \geq \overline{t}_i$. \end{itemize} \end{lemma} \begin{proof} Using the variation of constant formula, we can write \begin{equation*} E_i(t)=e^{-k_it}E_i(0)+e^{-k_it}\int_0^{t}e^{k_is}f_i(M(s))ds, \quad t\geq 0. \end{equation*} Then, we deduce that $$ 0 \leq E_{i}(t) \leq e^{-k_{i}t} E_{i}(0)+\frac{f_{i}(0)}{k_{i}}(1-e^{-k_{i}t}) \leq \max\big\{ E_{i}(0),\frac{f_{i}(0)}{k_{i}}\big\} . $$ Therefore, $E_i$ is bounded. Suppose that $E_{i}(t)>f_i(0)/k_i$. Then $$ E_i'(t)=-k_iE_i(t)+f_i(M(t))f_i(0)/k_i$. (i) Suppose that $E_i(t)>f_i(0)/k_i$, for all $t \geq 0$. Then, $E_i(t)$ is decreasing on $[0,+\infty)$, and $L_i:=\lim_{t \to +\infty}E_i(t)$ exists. Assume by contradiction that $L_i>f_i(0)/k_i$. Then $$ E_i'(t)+k_iE_i(t)=f_i(M(t)) \leq f_i(0), \quad t \geq 0. $$ By taking the limit in this last equation, we obtain $k_iL_i \leq f_i(0)$. This gives a contradiction. We conclude that $\lim_{t \to +\infty}E_i(t)=f_i(0)/k_i$. (ii) Suppose there exists $\overline{t}_i \geq 0$ such that $E_i(\overline{t}_i) = f_i(0)/k_i$. Then $$ E_i'(\overline{t}_i)=-f_i(0)+f_i(M(\overline{t}_i)) \leq 0. $$ Our objective is to prove that $E_i(t) \leq f_i(0)/k_i$, for all $t \geq \overline{t}_i$. If we suppose the existence of $\epsilon > 0$ such that $E_i(\overline{t}_i+\epsilon) > f_i(0)/k_i$, we obtain a contradiction, because the function $E_i(t)$ is strictly decreasing as long as $E_i(t) > f_i(0)/k_i$. \end{proof} Next, we state and prove a result regarding the boundedness property of the solutions of \eqref{sys.imp.suit}. In the rest of the paper, we suppose that $\delta >0$. The case $\delta=0$ should be treated separately, and so it will not be considered here. \begin{proposition}\label{P.4.3} Assume that \begin{equation}\label{Bound} \Big( 2\alpha \exp\Big( {-\tau\gamma\big(\frac{f_2(0)}{k_2}\big)}\Big) -1\Big)\beta\big(\frac{f_1(0)}{k_1}\big) < \delta. \end{equation} Then the solutions of system \eqref{sys.imp.suit} are bounded. \end{proposition} \begin{proof} A direct application of Lemma \ref{L.4.2} implies that $E_i$ are always bounded. Furthermore, it is not difficult to see that the boundedness of $N$ implies the boundedness of $M$. Then, we concentrate on the boundedness of $N$. By \eqref{Bound} and the continuity of $\gamma$ and $\beta$, we can take $\epsilon > 0$ small enough such that \begin{equation}\label{BoundEpsi} \Big( 2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) -1\Big)\beta\big(\frac{f_1(0)}{k_1}+\epsilon\big) < \delta. \end{equation} Lemma \ref{L.4.2} implies that there exits $\overline{t}_{\epsilon} \geq 0$ such that $E_i(t) \leq f_i(0)/k_i+\epsilon$, for all $t \geq \overline{t}_{\epsilon}$. Consider the function $$ Z_{\epsilon}(t)=N(t)+2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) \Big({ { \int_{t-\tau}^{t}{\beta(E_1(\theta))N(\theta)\, d\theta}}}\Big), \quad t \geq \overline{t}_{\epsilon}+\tau. $$ It follows that \begin{align*} Z_{\epsilon}'(t) &= N'(t)+2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) \Big(\beta(E_1(t))N(t) \\ & \quad -\beta(E_1(t-\tau))N(t-\tau)\Big),\\ &= - (\delta + \beta(E_1(t))) N( {t} ) \\ &\quad + 2\alpha \beta(E_1(t-\tau)) N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E_2(s))ds}\Big)\\ &\quad +2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) \beta(E_1(t))N(t)\\[0.5cm] &\quad -2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) \beta(E_1(t-\tau))N(t-\tau). \end{align*} This implies \begin{align*} &Z_{\epsilon}'(t)\\ &= -\Big[\delta-\Big(2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) -1\Big)\beta(E_1(t))\Big]N(t)\\ &\quad -2\alpha\Big(\exp({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}) -\exp(-{ \int_{t-\tau}^{t}{\gamma(E_2(s))\, ds}})\Big) \beta(E_1(t-\tau))N(t-\tau). \end{align*} As the function $\gamma$ is decreasing, we have $$ \exp\Big(-{ \int_{t-\tau}^{t}{\gamma(E_2(s))ds}}\Big) \leq \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big), \quad t \geq \overline{t}_{\epsilon}+\tau. $$ Then $$ Z_{\epsilon}'(t)\leq-\Big(\delta-\Big(2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big)-1\Big) \beta(E_1(t))\Big)N(t). $$ We have to consider two cases. Suppose that $$ 2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) \leq 1. $$ Then $ Z_{\epsilon}'(t) \leq 0$. Now, suppose that $$ 2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) > 1. $$ Since $\beta$ is increasing, we have $$ \beta(E_1(t))<\beta\big(\frac{f_1(0)}{k_1}+\epsilon\big). $$ Thanks to \eqref{BoundEpsi}, we obtain $$ \delta >\Big(2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2} +\epsilon\big)}\Big)-1\Big)\beta(E_1(t)). $$ This implies $ Z_{\epsilon}'(t) \leq 0$ for $t \geq \overline{t}_{\epsilon}+\tau$. We conclude that $Z_{\epsilon}$ is bounded. Consequently, $N$ is also bounded. \end{proof} \section{Existence of steady states}\label{S5} In this section, we study the existence of steady states of \eqref{sys.imp.suit}. Let $(\overline{N}, \overline{M}, \overline{E}_i)$ be a steady state of \eqref{sys.imp.suit}. Then, it satisfies \begin{equation}\label{steady_states} \begin{gathered} ( \delta + \beta(\overline{E}_1))\overline{N} = 2\alpha\beta(\overline{E}_1) \overline{N} e^{-\tau \gamma(\overline{E}_2)}, \\ \mu \overline{M} = K_N(\overline{E}_3)\overline{N} + 2K_P(\overline{E}_4)\beta(\overline{E}_1) \overline{N} e^{-\tau \gamma(\overline{E}_2)}, \\ k_i\overline{E}_i = f_i(\overline{M}). \end{gathered} \end{equation} One can easily see that $(0,0,f_i(0)/{k_i})$ is always a steady state (the trivial steady state). A nontrivial steady state $(\overline{N},\overline{M},\overline{E}_i) \neq (0,0,f_i(0)/{k_i})$, satisfies \begin{equation}\label{steady_states2} \begin{gathered} \delta =(2\alpha e^{-\tau\gamma(\overline{E}_2)}-1)\beta(\overline{E}_1), \\ \overline{N} = \frac{\mu\overline{M}}{K_{N}(\overline{E}_3)+2K_{P}(\overline{E}_{4})\beta(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)}}, \\ \overline{E}_{i} ={ \frac{f_{i}(\overline{M})}{k_{i}}.} \end{gathered} \end{equation} It is clear that, the existence and uniqueness of nontrivial steady state is equivalent to finding $\overline{M} > 0$, a solution of the equation \begin{equation}\label{SteadyStateM} \delta= \Big( 2\alpha \exp\Big( {-\tau\gamma\big(\frac{f_2(\overline{M})}{k_2}\big)}\Big) -1\Big)\beta(\frac{f_1(\overline{M})}{k_1}). \end{equation} \begin{proposition} Assume that \begin{equation}\label{PositiveSteady} \delta < \Big( 2\alpha \exp\Big( {-\tau\gamma\big(\frac{f_2(0)}{k_2}\big)}\Big) -1\Big) \beta\big(\frac{f_1(0)}{k_1}\big). \end{equation} Then there exists a unique nontrivial steady state $(\overline{N}, \overline{M}, \overline{E}_i)$ of \eqref{sys.imp.suit}. If \eqref{PositiveSteady} does not hold, then $(0,0,f_i(0)/{k_i})$ is the only steady state of \eqref{sys.imp.suit}. \end{proposition} \begin{proof} We define the function $$ \Psi(x)= \xi(x)\eta(x), \quad x \geq 0, $$ where $$ \xi(x)= 2\alpha \exp\Big( {-\tau\gamma(\frac{f_2(x)}{k_2})}\Big) -1, \quad \eta(x)=\beta\big(\frac{f_1(x)}{k_1}\big), \quad x \geq 0. $$ Then equation \eqref{SteadyStateM} becomes $\Psi(x)=\delta$ and inequality \eqref{PositiveSteady} can be written as $\delta<\Psi(0)$. Note that $\lim_{x \to +\infty}f_i(x)=0$ and $\beta(0)=0$. Then $\lim_{x \to +\infty}\Psi(x)=0$. We conclude that \eqref{SteadyStateM} has at least one solution if and only if $\delta <\Psi(0)$. To prove the uniqueness, we note that $\beta$ is increasing, $f_i$ and $\gamma$ are decreasing. Then, the functions $\xi$ and $\eta$ are decreasing and satisfy $\lim_{x \to +\infty}\xi(x)=2\alpha e^{-\tau \gamma(0)}-1$ and $\lim_{x \to +\infty}\eta(x)=0$. Firstly, we suppose that $2\alpha e^{-\tau \gamma(0)}-1\geq 0$. Then, $\xi(0)>0$. Consequently, $\Psi$ is positive, decreasing on $[0,+\infty)$ and satisfies $\lim_{x \to + \infty}\Psi(x)=0$. We conclude that \eqref{SteadyStateM} has a unique positive solution $\overline{M}$ if and only if $\delta < \Psi(0)$. Secondly, we suppose that $2\alpha e^{-\tau \gamma(0)}-1 <0$ and $\xi(0)>0$. Then, there exists a unique $\tilde{M} >0$ such that $\Psi(\tilde{M}) = 0$, $\Psi(x) > 0$ for $0\leq x < \tilde{M}$ and $\Psi(x) < 0$ for $x > \tilde{M}$. Consequently, $\Psi$ is positive and decreasing on the interval $[0,\tilde{M}]$ with $\Psi(\tilde{M})=0$. Then \eqref{SteadyStateM} has a unique solution $\overline{M} \in (0,\tilde{M})$ if and only if $\delta < \Psi(0)$. The existence and uniqueness of $\overline{E}_i$ and $\overline{N}$ follow immediately from \eqref{steady_states2}. This completes the proof. \end{proof} Note that condition \eqref{PositiveSteady} is equivalent to \begin{equation}\label{ineglty2.imp.suit} \begin{gathered} 1 \geq \alpha>\alpha_{\rm min} := \frac{\delta+\beta(f_1(0)/k_1)}{2\beta(f_1(0)/k_1)} = \frac{1}{2}+ \frac{\delta}{2\beta(f_1(0)/k_1)},\\ 0\leq\tau<\tau_{\rm max} := \frac{1}{\gamma(f_2(0)/k_2)} \ln\Big(\frac{2\alpha\beta(f_1(0)/k_1)}{\delta+\beta(f_1(0)/k_1)}\Big). \end{gathered} \end{equation} In the next section, we analyze the asymptotic behavior of the solutions of system \eqref{sys.imp.suit} by studying the asymptotic stability of its steady states. \section{Global asymptotic stability of trivial steady state}\label{S6} We assume, throughout this section, that the function $\beta$ is continuously differentiable on $[0, +\infty)$. We begin by establishing the global asymptotic stability of the trivial steady state $(0,0,f_{i}(0)/k_{i})$. First let us recall a useful lemma (see Gopalsamy \cite{Gopalsamy1992}), that will allow us to establish the next result. \begin{lemma}\label{Lemma 6.1} Let $f:(a,+\infty) \to \mathbb{R}$, $a \in \mathbb{R}$, be a differentiable function. If $\lim_{t \to +\infty}f(t)$ exists and $f'(t)$ is uniformly continuous on $(a,+\infty)$, then \[ \lim_{t \to +\infty}f'(t)=0. \] \end{lemma} \begin{lemma}\label{lemme6.2} Let $(N(t),M(t),E_{i}(t))$ be a bounded solution of \eqref{sys.imp.suit}. Then, the following three statements are equivalent $$ \lim_{t \to +\infty}N(t)=0, \quad \lim_{t \to +\infty}M(t)=0,\quad \lim_{t \to +\infty}E_i(t)=f_i(0)/k_i, \; i=1,2,3,4. $$ \end{lemma} \begin{proof} We begin by proving that $\lim_{t \to +\infty}N(t)=0$ if and only if $\lim_{t \to +\infty}M(t)=0$. We first assume that $\lim_{t \to +\infty}N(t)=0$. We have $M'(t)=-\mu M(t)+F(t)$, with \begin{align*} F(t)&= K_N(E_3(t))N(t) + 2K_P(E_4(t))\beta(E_1(t-\tau))\\ &\quad\times N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E_2(s))\,ds}\Big). \end{align*} Then $\lim_{t \to +\infty}F(t)=0$. Using the variation of constant formula, we can write \[ M(t)=e^{-\mu t}M(0)+e^{-\mu t}\int_0^{t}e^{\mu s}F(s)ds, \quad t\geq 0. \] Let $\epsilon >0$ be fixed. Since $N(t)$ tends to zero when $t$ tends to $+\infty$, there exits $t_\epsilon >0$ such that \[ F(t)<\frac{\mu \epsilon}{2}, \quad e^{-\mu t}\Big(M(0)+\int_0^{t_\epsilon}e^{\mu s}F(s)ds\Big) < \frac{\epsilon}{2}, \quad \text{for } t\geq t_\epsilon. \] Then, for $t \geq t_\epsilon$, we have \[ M(t)\leq e^{-\mu t}\Big(M(0)+\int_0^{t_\epsilon}e^{\mu s}F(s)ds\Big) +e^{-\mu t}\Big(\int_{t_\epsilon}^{t}e^{\mu s}F(s)ds\Big), \] with \[ e^{-\mu t}\Big(\int_{t_\epsilon}^{t}e^{\mu s}F(s)ds\Big) \leq \frac{\epsilon}{2} \big( 1-e^{\mu(t_\epsilon-t)}\big) \leq \frac{\epsilon}{2}. \] Consequently, $M(t)<\epsilon$ for $t\geq t_\epsilon$. We have proved that $\lim_{t \to +\infty}M(t)=0$. Secondly, we assume that $\lim_{t \to +\infty}M(t)=0$. Then the solution $(N, M, E_i)$ is bounded, and the derivative $(N', M', E'_i)$ is also bounded. Furthermore, $M'(t)$ is differentiable for $t>\tau$, and since $N$, $M$, $E_i$, $N'$, $M'$, $E'_i$ are bounded for $t>\tau$, $M''$ is bounded. Then $M'$ is uniformly continuous. Consequently, Lemma \ref{Lemma 6.1} implies that $\lim_{t \to +\infty}M'(t)=0$. From the equation satisfied by $M$, we deduce that $\lim_{t \to +\infty}F(t)=0$. In particular, $\lim_{t \to +\infty}K_N(E_3(t))N(t)=0$. Suppose that $\lim_{t \to +\infty}K_N(E_3(t))=0$. That means that $\lim_{t \to +\infty}E_3(t)=0$. Then, from the equation satisfied by $E_3$, we deduce that $\lim_{t \to +\infty}f_3(M(t))=0$. This gives a contradiction because $\lim_{t \to +\infty}M(t)=0$ and $f_3(0)>0$. Then, we conclude that $\lim_{t \to +\infty}N(t)=0$.\\ Thirdly, we assume that $\lim_{t \to +\infty}M(t)=0$ and we will prove that $\lim_{t \to +\infty}E_i(t)=f_i(0)/k_i$. Thanks to Lemma \ref{L.4.2}, it suffices to prove the result for the case $E_i(t) \bar{t}$. Without loss of generality, we can choose $\bar{t}=0$. We put $F_i(t)=f_i(0)/k_i-E_i(t)$ and $G_i(t)=f_i(0)-f_i(M(t))$. Then, $F_i$ satisfies the following differential equation \begin{equation*} F'_i(t)=-k_iF_i(t)+ G_i(t), \end{equation*} with $F_i(t)>0$ and $G_i(t)>0$ for all $t \geq 0$ and $\lim_{t \to +\infty}G_i(t)=0$. Then, using the same argument as in the first part of this proof, we obtain $\lim_{t \to +\infty}F_i(t)=0$. This means that $\lim_{t \to +\infty}E_i(t)=f_i(0)/k_i$. Finally, we suppose that $\lim_{t \to +\infty}E_i(t)=f_i(0)/k_i$. Then by Lemma \ref{Lemma 6.1}, we have $\lim_{t \to +\infty}E'_i(t)=0$. We deduce that $\lim_{t \to +\infty}f_i(M(t))=f_i(0)$. As the function $f_i$ is continuous and strictly decreasing from $[0,+\infty)$ into $(0,f_i(0)]$, we conclude that $\lim_{t \to +\infty}M(t)=0$. This completes the proof. \end{proof} Next, we prove a result dealing with the global asymptotic stability of the trivial steady state $(0,0,f_{i}(0)/k_{i})$. \begin{theorem}\label{globalstabTrivialSteady} Assume that \eqref{Bound} holds. That is to say that $(0,0,f_{i}(0)/k_{i})$ is the only steady state. Then, all solutions $(N(t),M(t),E_{i}(t))$ of \eqref{sys.imp.suit} converge to $(0,0,f_{i}(0)/k_{i})$, $i=1,2,3,4$. \end{theorem} \begin{proof} As in the proof of Proposition \ref{P.4.3}, we take $\epsilon > 0$ small enough such that \begin{equation*} \Big( 2\alpha \exp\Big({-\tau \gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)}\Big) -1\Big) \beta\big(\frac{f_1(0)}{k_1}+\epsilon\big) < \delta, \end{equation*} and $\overline{t}_{\epsilon} \geq 0$ such that $E_i(t) \leq f_i(0)/k_i+\epsilon$, for all $t \geq \overline{t}_{\epsilon}$. Consider the functional $ V_{\epsilon} : (C([\overline{t}_{\epsilon},\overline{t}_{\epsilon}+\tau], \mathbb{R}_{+}))^{6} \to \mathbb{R}_{+}$, Defined by \begin{gather*} \varPhi =(\varphi,\psi,\chi_1,\chi_2,\chi_3,\chi_{4}),\\ \begin{aligned} V_{\epsilon}(\varPhi)&= \varphi(\overline{t}_{\epsilon}+\tau) +2\alpha \exp\Big(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)\Big)\\ &\quad\times \int_{-\tau}^{0}\beta(\chi_1(\theta+\overline{t}_{\epsilon}+\tau)) \varphi(\theta+\overline{t}_{\epsilon}+\tau)\, d\theta. \end{aligned} \end{gather*} The composition with the solution $X(t):=(N(t),M(t),E_i(t) )$ of equation \eqref{sys.imp.suit} leads, for $t \geq \overline{t}_{\epsilon}+\tau$, to the function \begin{equation*} t\mapsto V_{\epsilon}(X_{t}) =N(t)+2\alpha \exp\Big(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)\Big) \int_{t-\tau}^{t}\beta(E_1(s))N(s)ds. \end{equation*} Then, the derivative along the solution of system \eqref{sys.imp.suit} gives \begin{align*} &\frac{d}{dt}V_{\epsilon}(X_{t})\\ &= N'(t) +2\alpha \exp\Big(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)\Big) \big[\beta(E_1(t))N(t)-\beta(E_1(t-\tau))N(t-\tau)\big],\\ &= -(\delta+\beta(E_1(t)))N(t)\\ &\quad +2\alpha\beta(E_1(t-\tau))N(t-\tau) \exp\Big(-{ \int_{t-\tau}^{t}{\gamma(E_2(s))\, ds}}\Big)\\ &\quad +2\alpha \exp(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)) [\beta(E_1(t))N(t)-\beta(E_1(t-\tau))N(t-\tau)],\\ &= -\Big(\delta-\big[2\alpha \exp(-\tau\gamma \big(\frac{f_2(0)}{k_2}+\epsilon\big))-1\big]\beta(E_1(t)))N(t)\\ &\quad -2\alpha\big[\exp(-\tau\gamma\big(\frac{f_2(0)}{k_2} +\epsilon\big))-\exp(-{ \int_{-\tau}^{0}{\gamma(E_2(t+s))\, ds}})\big]\\ &\quad \times\beta(E_1(t-\tau))N(t-\tau). \end{align*} Let $-\tau \leq s \leq 0$. Since $t \geq \overline{t}_{\epsilon}+\tau$, then $E_i(t+s) < { \frac{f_{i}(0)}{k_{i}}}+\epsilon$. Consequently, $-\gamma(E_2(s+t))<-\gamma({ \frac{f_2(0)}{k_2}}+\epsilon)$ and $\beta(E_1(t))<\beta({ \frac{f_1(0)}{k_1}} + \epsilon)$. This implies \begin{gather*} \exp\Big(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)\Big) > \exp\Big(-{ \int_{-\tau}^{0}{\gamma(E_2(s+t))\, ds}}\Big), \\ \delta > \Big[2\alpha \exp\Big(-\tau\gamma\big(\frac{f_2(0)}{k_2}+\epsilon\big)\Big) -1\Big]\beta(E_1(t)). \end{gather*} Thus, $$ \dot{V}_\epsilon(\varPhi)\leq 0, \quad \text{for all } \varPhi \in (C([\overline{t}_{\epsilon}, \overline{t}_{\epsilon}+\tau],\mathbb{R}_{+}))^{6}, $$ where $\dot{V}_\epsilon$ is the derivative of $V_\epsilon$ along the solutions of \eqref{sys.imp.suit}. Now, let $$ S:=\{\varPhi \in (C([\overline{t}_{\epsilon},\overline{t}_{\epsilon} +\tau],\mathbb{R}_{+}))^{6}: \dot{V}_\epsilon(\varPhi)=0\}. $$ We deduce that $$ S=\{\varPhi \in (C([\overline{t}_{\epsilon},\overline{t}_{\epsilon}+\tau], \mathbb{R}_{+}))^{6}: \varphi(\overline{t}_{\epsilon}+\tau) =\varphi(\overline{t}_{\epsilon})=0\}. $$ We also consider the set $\Omega$, defined as the largest set in $S$ which is invariant with respect to system \eqref{sys.imp.suit}. Let $X_t$ be a solution of \eqref{sys.imp.suit} associated with an initial condition $\varPhi \in \Omega$. Then, $X_t \in \Omega$ for all $t \geq \overline{t}_{\epsilon}+\tau$ is equivalent to $N(t)=0$ for all $t \geq \overline{t}_{\epsilon}+\tau$. Consequently, $$ \Omega=\{0\}\times(C([\overline{t}_{\epsilon}, \overline{t}_{\epsilon}+\tau],\mathbb{R}_{+}))^{5}. $$ From Hale and Verduyn Lunel \cite[page 143]{Hale1993S}, all bounded solutions $X_t$ of \eqref{sys.imp.suit} converge to $\Omega$ as $t$ tends to $+\infty$. From Proposition \ref{Bound}, all solutions of \eqref{sys.imp.suit} are bounded provided that \eqref{Bound} holds. Then, all solutions $X_t$ converge to $\Omega$. We deduce that for all $\varPhi \in (C([\overline{t}_{\epsilon},\overline{t}_{\epsilon}+\tau], \mathbb{R}_{+}))^{6}$, $$ \lim_{t \to +\infty}N(t)=0. $$ Then, from Lemma \ref{lemme6.2}, we conclude that \[ \lim_{t \to +\infty}M(t)=0 \quad \text{and} \quad \lim_{t \to +\infty}E_i(t)=f_i(0)/k_i, \quad i=1,2,3,4.\] This completes the proof. \end{proof} Next, we linearize system \eqref{sys.imp.suit} about its steady states, and we determine the characteristic equation. Let $(\overline{N},\overline{M},\overline{E}_i)$ be a steady state of \eqref{sys.imp.suit}. We set $$ X(t) = N(t)-\overline{N},\quad Y(t)=M(t)-\overline{M},\quad Z_i(t)=E_i(t)-\overline{E}_i. $$ The linearized system of \eqref{sys.imp.suit} around $(\overline{N},\overline{M},\overline{E}_i)$ is \begin{equation}\label{sys.liner.gen} \begin{aligned} X'(t) & =-(\delta+\beta(\overline{E}_1))X(t)-\beta'(\overline{E}_1)\overline{N}Z_1(t)\\ &\quad +2\alpha\beta(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)}X(t-\tau) +2\alpha\overline{N}\beta'(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)} Z_1(t-\tau)\\ & \quad -2\tau\alpha\overline{N}\beta(\overline{E}_1)\gamma'(\overline{E}_2) e^{-\tau\gamma(\overline{E}_2)}{ \int_{-\tau}^{0}Z_2(t+s)\, ds},\\ Y'(t) & =-\mu Y(t)+K_{N}(\overline{E}_3)X(t)+K_{N}'(\overline{E}_3) \overline{N}Z_3(t)\\ & \quad +2K_{P}(\overline{E}_{4})\beta(\overline{E}_1) e^{-\tau\gamma(\overline{E}_2)}X(t-\tau) +2K_{P}'(\overline{E}_{4})\beta(\overline{E}_1) \overline{N}e^{-\tau\gamma(\overline{E}_2)}Z_{4}(t)\\ & \quad +2K_{P}(\overline{E}_{4})\beta'(\overline{E}_1)\overline{N} e^{-\tau\gamma(\overline{E}_2)}Z_1(t-\tau)\\ & \quad -2\tau K_{P}(\overline{E}_{4})\beta(\overline{E}_1) \overline{N}\gamma'(\overline{E}_2) e^{-\tau\gamma(\overline{E}_2)}{ \int_{-\tau}^{0}Z_2(t+s)\, ds},\\ Z_{i}'(t) & =-k_{i}Z_{i}(t)+f_{i}'(\overline{M})Y(t). \end{aligned} \end{equation} The above system has the form \begin{equation}\label{sys.liner.gen.matrices} U'(t)=AU(t)+BU(t-\tau)+C({ \int_{-\tau}^{0}U(t+s)\, ds}), \end{equation} with $U(t)=(X(t),Y(t),Z_{i}(t))^{T}\in\mathbb{R}^{6}$, where \begin{gather*} A=\begin{pmatrix} -(\delta+\beta(\overline{E}_1)) & 0 & -\beta'(\overline{E}_1)\overline{N} & 0 & 0 & 0\\ K_{N}(\overline{E}_3) & -\mu & 0 & 0 & K_{N}'(\overline{E}_3)\overline{N} & \overline{N} \overline{H}_4\\ 0 & f_1'(\overline{M}) & -k_1 & 0 & 0 & 0\\ 0 & f_2'(\overline{M}) & 0 & -k_2 & 0 & 0\\ 0 & f_3'(\overline{M}) & 0 & 0 & -k_3 & 0\\ 0 & f_{4}'(\overline{M}) & 0 & 0 & 0 & -k_{4} \end{pmatrix}, \\ B=\begin{pmatrix} H_0 & 0 & \overline{N}\,H_1 & 0 & 0 & 0\\ \overline{H}_0 & 0 & \overline{N}\,\overline{H}_1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} ,\quad C=\begin{pmatrix}0 & 0 & 0 & \overline{N}\,H_2 & 0 & 0\\ 0 & 0 & 0 & \overline{N}\,\overline{H}_2 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \end{gather*} with \begin{gather*} H_0 := 2\alpha \beta(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)} \geq 0, \\ H_1 := \frac{dH_0}{d\overline{E}_1} = 2\alpha \beta'(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)} \geq 0, \\ H_2 := \frac{dH_0}{d\overline{E}_2} = -2\tau \alpha \beta(\overline{E}_1)\gamma' (\overline{E}_2) e^{-\tau\gamma(\overline{E}_2)} \geq 0, \end{gather*} and \begin{gather*} \overline{H}_0 := 2K_{P}(\overline{E}_{4})\beta(\overline{E}_1) e^{-\tau\gamma(\overline{E}_2)} \geq 0, \\ \overline{H}_1 := \frac{d\overline{H}_0}{d\overline{E}_1} = 2K_{P}(\overline{E}_{4})\beta'(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)} \geq 0, \\ \overline{H}_2 := \frac{d\overline{H}_0}{d\overline{E}_2} = -2\tau K_{P}(\overline{E}_{4})\beta(\overline{E}_1)\gamma' (\overline{E}_2) e^{-\tau\gamma(\overline{E}_2)} \geq 0,\\ \overline{H}_4 := \frac{d\overline{H}_0}{d\overline{E}_{4}} = 2K'_{P}(\overline{E}_{4})\beta(\overline{E}_1)e^{-\tau\gamma(\overline{E}_2)} \geq 0. \end{gather*} The relationship between the expressions $H_0$ and $\overline{H}_0$ is given by \[ \alpha \overline{H}_0=K_{P}(\overline{E}_{4})H_0. \] The characteristic equation associated to the steady state $(\overline{N},\overline{M},\overline{E}_i)$ is \begin{equation}\label{thefirst.charac.equa} \Delta(\lambda)=\det\Big(\lambda I-A-Be^{-\lambda\tau}-{ C\int_{-\tau}^{0} e^{\lambda\theta}\, ds}\Big)=0. \end{equation} Then, we have the following result. \begin{theorem} The trivial steady state of \eqref{sys.imp.suit} is unstable when \eqref{PositiveSteady} holds. \end{theorem} \begin{proof} When $(\overline{N},\overline{M},\overline{E}_i) = (0,0,f_{i}(0)/k_{i})$ system \eqref{sys.liner.gen.matrices} becomes \begin{equation}\label{sys.liner.matrices2} U'(t)=AU(t)+BU(t-\tau), \end{equation} with $U(t)=(X(t),Y(t),Z_i(t))^{T}\in\mathbb{R}^{6}$, \[ A=\begin{pmatrix}-(\delta+\beta(f_1(0)/k_1)) & 0 & 0 & 0 & 0 & 0\\ K_{N}( f_3(0)/k_3 ) & -\mu & 0 & 0 & 0 & 0\\ 0 & f_1'(0) & -k_1 & 0 & 0 & 0\\ 0 & f_2'(0) & 0 & -k_2 & 0 & 0\\ 0 & f_3'(0) & 0 & 0 & -k_3 & 0\\ 0 & f_{4}'(0) & 0 & 0 & 0 & -k_{4} \end{pmatrix}, \] and \[ B=2\beta( \frac{f_1(0)}{k_1}) \exp\Big( {-\tau\gamma\big(\frac{f_2(0)}{k_2}\big)}\Big) \begin{pmatrix} \alpha & 0 & 0 & 0 & 0 & 0\\ K_{P}(f_{4}(0)/{k_{4}}) & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}. \] The characteristic equation \eqref{thefirst.charac.equa} becomes \begin{equation*} \Delta(\lambda)=\det(\lambda I-A-Be^{-\lambda\tau})=0. \end{equation*} Then \[ \Delta(\lambda)=(\lambda+\mu) \prod_{i=1}^4 (\lambda+k_{i})\overline{\Delta}(\lambda), \] where \[ \overline{\Delta}(\lambda)=\lambda+\delta+\beta({ \frac{f_1(0)}{k_1}}) -2\alpha\beta({ \frac{f_1(0)}{k_1}}) \exp\Big(-\tau(\lambda+\gamma({ \frac{f_2(0)}{k_2}}))\Big). \] The eigenvalues of \eqref{sys.liner.matrices2} are $\lambda=-\mu<0$, $\lambda=-k_{i}<0,\;i\in{1,2,3,4}$ and roots of the equation $\overline{\Delta}(\lambda)=0$. Let $\lambda \in \mathbb{R}$. We have \begin{gather*} \overline{\Delta}'(\lambda)=1+2\alpha\tau\beta({ \frac{f_1(0)}{k_1}}) \exp\Big(-\tau(\lambda+\gamma({ \frac{f_2(0)}{k_2}}))\Big)>0, \\ \overline{\Delta}(0)=\delta-2\alpha\Big( \exp\Big(-\tau\gamma({ \frac{f_2(0)}{k_2}})\Big)-1\Big) \beta\big({ \frac{f_1(0)}{k_1}}\big), \\ \lim_{\lambda\to+\infty} \overline{\Delta}(\lambda)=+\infty. \end{gather*} Thanks to \eqref{PositiveSteady}, we have $\overline{\Delta}(0) < 0$. Thus, there exists $\lambda_0>0$ such that $\overline{\Delta}(\lambda_0)=0$. Hence, the instability of the trivial steady state holds. \end{proof} The last theorem completes the global asymptotic stability of the trivial steady state $(0,0,f_{i}(0)/k_{i})$ obtained in Theorem \ref{globalstabTrivialSteady}, and allows us to entirely determine its dynamics. \section{Local asymptotic stability of the positive steady state}\label{S7} We assume throughout this section, that condition \eqref{PositiveSteady} is satisfied, or equivalently \eqref{ineglty2.imp.suit}, to ensure the existence and uniqueness of the positive steady state $(\overline{N},\overline{M},\overline{E}_i)$ of \eqref{sys.imp.suit}. The nature of the characteristic equation associated to the linearized system around $(\overline{N},\overline{M},\overline{E}_i)$ induces some technical difficulties. To avoid these difficulties, we make the following assumption \[ k_i=k, \quad f_i=f, \quad E_i=E, \quad i=1,2,3,4. \] Then system \eqref{sys.imp.suit} becomes \begin{equation}\label{RedDela} \begin{aligned} N'(t) &= - (\delta + \beta(E(t))) N( {t} ) \\ &\quad + 2\alpha\beta(E(t-\tau)) N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E(s))\,ds}\Big),\\ M'(t) &= - \mu M( {t} ) + K_N(E(t))N(t) \\ &\quad + 2K_P(E(t))\beta(E(t-\tau)) N(t-\tau) \exp \Big( - \int_{t-\tau}^{t} { \gamma(E(s))\,ds}\Big),\\ E'(t) &= -kE(t) + f(M(t)). \end{aligned} \end{equation} The linearized system of \eqref{RedDela} around the positive steady state $(\overline{N},\overline{M},\overline{E})$ has the form \begin{equation}\label{eq.sys.lin} U'(t)=AU(t)+BU(t-\tau)+C\Big({ \int_{-\tau}^{0}U(t+s)\, ds}\Big), \end{equation} with $U(t)=(X(t),Y(t),Z(t))^{T}\in\mathbb{R}^{3}$, \begin{gather*} A=\begin{pmatrix}-(\delta+\beta(\overline{E})) & 0 & -\beta'(\overline{E})\overline{N} \\ K_{N}(\overline{E}) & -\mu & K_{N}'(\overline{E})\overline{N}+2K_{P}'(\overline{E})\beta(\overline{E})\overline{N}e^{-\tau\gamma(\overline{E})} \\ 0 & f'(\overline{M}) & -k \end{pmatrix},\\ B=\begin{pmatrix}2\alpha\beta(\overline{E})e^{-\tau\gamma(\overline{E})} & 0 & 2\alpha\overline{N}\beta'(\overline{E})e^{-\tau\gamma(\overline{E})}\\ 2K_{P}(\overline{E})\beta(\overline{E})e^{-\tau\gamma(\overline{E})} & 0 & 2K_{P}(\overline{E})\beta'(\overline{E})\overline{N}e^{-\tau\gamma(\overline{E})}\\ 0 & 0 & 0 \end{pmatrix},\\ C=\begin{pmatrix}0 & 0 & -2\tau\alpha\overline{N}\beta(\overline{E}) \gamma'(\overline{E})e^{-\tau\gamma(\overline{E})}\\ 0 & 0 & -2\tau K_{P}(\overline{E})\beta(\overline{E})\overline{N} \gamma'(\overline{E})e^{-\tau\gamma(\overline{E})}\\ 0 & 0 & 0 \end{pmatrix}. \end{gather*} The associated characteristic equation becomes \begin{equation}\label{general.charac.equa} \Delta(\lambda) =\det\Big(\lambda I-A-Be^{-\lambda\tau}-{ C\int_{-\tau}^{0} e^{\lambda\theta}\, ds}\Big)=0. \end{equation} The condition \eqref{ineglty2.imp.suit} becomes \begin{equation}\label{ineglty3} \begin{gathered} 1 \geq \alpha>\alpha_{\rm min}:=\frac{1}{2}+ \frac{\delta}{2\beta(f(0)/k)},\\ 0\leq\tau<\tau_{\rm max}:= \frac{1}{\gamma(f(0)/k)} \ln\Big(\frac{2\alpha\beta(f(0)/k)}{\delta+\beta(f(0)/k)}\Big). \end{gathered} \end{equation} First, let us suppose that $\tau=0$. Then \eqref{eq.sys.lin} becomes \begin{equation}\label{eq.sys.lin.2} U'(t)=(A+B)U(t), \end{equation} where $$ A+B=\begin{pmatrix}-(\delta-(2\alpha-1)\beta(\overline{E})) & 0 & (2\alpha-1)\beta'(\overline{E})\overline{N}\\ \Lambda(\overline{E}) & -\mu & \Lambda'(\overline{E})\overline{N}\\ 0 & f'(\overline{M}) & -k \end{pmatrix}, $$ with $\Lambda(\overline{E})=K_{N}(\overline{E})+2K_{P}(\overline{E}) \beta(\overline{E})$. In fact, under condition \eqref{ineglty3}, the steady state $(\overline{N},\overline{M},\overline{E})$ satisfies \[ (2\alpha -1)\beta(\overline{E})=\delta, \quad \Lambda(\overline{E})\overline{N}=\mu \overline{M},\quad k\overline{E}=f(\overline{M}). \] Then the characteristic equation of \eqref{eq.sys.lin.2}, \begin{equation*} \det(\lambda I-A-B)=0, \end{equation*} becomes $$ \lambda\left[(\lambda+\mu)(\lambda+k)-f'(\overline{M}) \Lambda'(\overline{E})\overline{N}\right]-(2\alpha-1) \Lambda(\overline{E})f'(\overline{M})\beta'(\overline{E})\overline{N}=0. $$ This is equivalent to \begin{equation}\label{eq.charac.param.routh} \lambda^{3}+a_1\lambda^{2}+a_2\lambda+a_3=0, \end{equation} where \begin{gather*} a_1=\mu+k >0, \quad a_2=\mu k-f'(\overline{M})\Lambda'(\overline{E})\overline{N}>0,\\ a_3=-(2\alpha-1)\Lambda(\overline{E})f'(\overline{M})\beta'(\overline{E}) \overline{N}>0. \end{gather*} We have \[ a_1a_2-a_3=\mu k(\mu+k)-\overline{N}f'(\overline{M}) [(\mu+k)\Lambda'(\overline{E}) -(2\alpha-1)\Lambda(\overline{E}) \beta'(\overline{E})]. \] Then, by applying the Ruth-Hurwitz criterion, we obtain the following lemma. \begin{lemma}\label{lemma7.1} All roots of \eqref{eq.charac.param.routh} have negative real parts if and only if \begin{equation}\label{criterion.routh} \mu k(\mu+k)-\overline{N}f'(\overline{M})\left[(\mu+k)\Lambda'(\overline{E}) -(2\alpha-1)\Lambda(\overline{E})\beta'(\overline{E})\right]>0. \end{equation} \end{lemma} We also have the following lemma. \begin{lemma} Assume that \begin{equation}\label{NewCond} \mu+k>\frac{(2\alpha-1)\Lambda(\beta^{-1}(\delta/{2\alpha-1} )) \beta'(\beta^{-1}(\delta/{2\alpha-1} ))}{\Lambda'(\beta^{-1}(\delta/{2\alpha-1} ))}, \end{equation} Then, all roots of \eqref{eq.charac.param.routh} have negative real parts. \end{lemma} \begin{proof} Since $f'(\overline{M})<0$, $\Lambda'(\overline{E})>0$ and $\beta'(\overline{E})>0$, the hypothesis \[ (\mu+k)\Lambda'(\overline{E}) >(2\alpha-1)\Lambda(\overline{E})\beta'(\overline{E}) \] implies that \eqref{criterion.routh} is satisfied. Furthermore, we have $\overline{E}=\beta^{-1}(\delta/{2\alpha-1} )$. Then, \eqref{NewCond} implies \eqref{criterion.routh}. Consequently, if \eqref{NewCond} is satisfied, then all roots of \eqref{eq.charac.param.routh} have negative real parts. \end{proof} \begin{theorem} Assume that \eqref{PositiveSteady} and \eqref{criterion.routh} hold. Then there exists $\tau^* \in [0,\tau_{\rm max})$ such that the positive steady state $(\overline{N},\overline{M},\overline{E})$ is locally asymptotically stable for $\tau \in [0,\tau^*)$. \end{theorem} \begin{proof} A direct application of Lemma \ref{lemma7.1} implies that $(\overline{N},\overline{M},\overline{E})$ is locally asymptotically stable when $\tau=0$. Furthermore, $\Delta:=\Delta(\lambda, \tau)$ given by \eqref{general.charac.equa} is analytic in $\lambda$ and $\tau$. Then, as $\tau$ varies the zeros of $\lambda \mapsto \Delta(\lambda,\tau)$ stay in the open left half-plane for $\tau$ small enough. If instability occurs for a particular value of $\tau$, a characteristic root must intersect the imaginary axis. Then, there exists $\tau^* \in [0,\tau_{\rm max})$ such that for $\tau \in [0,\tau^*)$, all the roots of \eqref{general.charac.equa} have negative real parts. \end{proof} \section{Numerical illustrations}\label{S8} In this section, we perform some numerical simulations to illustrate the behavior of the steady states. Let us choose the functions $\beta$, $\gamma$ , $K_P$, $K_N$ and the negative feedback $f$ as follows \begin{gather*} \beta (E) = \frac{\hat{\beta} E}{1+E}, \quad \gamma(E)=\frac{\gamma_0}{1+E^{a}}, \quad K_P(E)=\frac{\hat{K}_{P}E}{1+E}, \\ K_N(E)=\frac{\hat{K}_{N}E}{1+E}, \quad f(M)=\frac{f_0}{1+M^b}, \end{gather*} with \[ \delta + \hat{\beta} +\hat{K}_{N} \leq 1 \quad\text{and}\quad \alpha + \hat{K}_{P} \leq 1. \] From \eqref{ineglty2.imp.suit} the positive steady state $(\overline{N}, \overline{M}, \overline{E})$ exists if and only if \begin{gather*} 1 \geq \alpha>\alpha_{\rm min}:= \frac{1}{2}+ \frac{\delta}{2\beta(f_1(0)/k_1)},\\ 0\leq\tau<\tau_{\rm max}:= \frac{1}{\gamma(f_2(0)/k_2)} \ln\Big(\frac{2\alpha\beta(f_1(0)/k_1)}{\delta+\beta(f_1(0)/k_1)}\Big). \end{gather*} We fix all the parameters except the delay $\tau$. The values of the fixed parameters are $\delta=0.08$, $\mu=0.05$, $k=0.6$, $\alpha=0.8$, $\hat{\beta}=0.8$, $\gamma_0=0.2$, $a=3$, $\hat{K}_{P}= 0.18$, $\hat{K}_{N}= 0.1$, $f_0=1$, $b=7$. With these values, we have $\tau_{\rm max}=9.052$ and $\alpha_{\rm min}= 0.58$. Then, the above condition of existence of positive steady state becomes $0 \leq \tau < 9.052$. \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig2} % triv_gas1.eps \end{center} \caption{Global asymptotic stability of the trivial steady state $(\overline{N}=0,\overline{M}=0)$. When $\tau =9.5 > \tau_{\rm max}=9.052$, the trivial steady state is the only steady state and it is globally asymptotically stable.} \label{Fig2} \end{figure} \begin{figure}[htb] \begin{center} \includegraphics[width=0.6\textwidth]{fig3} %non_triv_las1_02.eps \end{center} \caption{Local asymptotic stability of the positive steady state $(\overline{N},\overline{M})$. When $\tau=1 < \tau_{\rm max}=9.052$, the positive steady state is locally asymptotically stable} \label{Fig3} \end{figure} \begin{thebibliography}{10} \bibitem{AdamsonAJM1996} John~W. Adamson; \emph{Regulation of red blood cell production}, The American Journal of Medicine \textbf{101} (1996), no.~2, 4S -- 6S. \bibitem{AdimyCrausteMMNP2012} M.~Adimy, F.~Crauste; \emph{Delay differential equations and autonomous oscillations in hematopoietic stem cell dynamics modeling}, Mathematical Modelling of Natural Phenomena \textbf{7} (2012), 1--22. \bibitem{AdimyCrausteDCDSB2007} Mostafa Adimy, Fabien Crauste; \emph{Modeling and asymptotic stability of a growth factor-dependent stem cell dynamics model with distributed delay}, Discrete Contin. Dyn. Syst. Ser. B \textbf{8} (2007), no.~1, 19--38. \bibitem{AdimyCrausteMCM2009} Mostafa Adimy, Fabien Crauste; \emph{Mathematical model of hematopoiesis dynamics with growth factor-dependent apoptosis and proliferation regulations}, Math. Comput. Modelling \textbf{49} (2009), no.~11-12, 2128--2137. \bibitem{AdimyCrausteRuanSIAP2005} Mostafa Adimy, Fabien Crauste, Shigui Ruan; \emph{A mathematical study of the hematopoiesis process with applications to chronic myelogenous leukemia}, SIAM J. Appl. Math. \textbf{65} (2005), no.~4, 1328--1352. \bibitem{AdimyCrausteRuanBMB2006} Mostafa Adimy, Fabien Crauste, Shigui Ruan; \emph{Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases}, Bull. Math. Biol. \textbf{68} (2006), no.~8, 2321--2351. \MR{2293845 (2007k:92034)} \bibitem{BelairMackeyMahaffyMB1995} Jacques {B\'elair}, Michael~C. {Mackey}, Joseph~M. {Mahaffy}; \emph{Age-structured and two-delay models for erythropoiesis.}, Math. Biosci. \textbf{128} (1995), no.~1-2, 317--346 (English). \bibitem{BurnsTannock1970} Fredric J. Burns, Ian F. Tannock; \emph{On existence of a g$_0$-phase in cell cycle}, Cell and Tissue Kinetics \textbf{3} (1970), no.~4, 321--334. \bibitem{Gopalsamy1992} K. Gopalsamy; \emph{Stability and oscillations in delay differential equations of population dynamics}, Mathematics and its Applications, vol.~74, Kluwer Academic Publishers Group, Dordrecht, 1992. \bibitem{Hale1993S} Jack K. Hale, Sjoerd M. Verduyn Lunel; \emph{Introduction to functional differential equations}, Applied Mathematical Sciences 99, Springer-Verlag, New York, 1993 (en). \bibitem{HaurieDaleMackeyBLOOD1998} Caroline Haurie, David C. Dale, Michael C. Mackey; \emph{Cyclical neutropenia and other periodic hematological disorders: A review of mechanisms and mathematical models}, Blood \textbf{{92}} ({1998}), no.~{8}, {2629--2640}. \bibitem{KouryBondurant1990} Mark J. Koury, Maurice C. Bondurant; \emph{Erythropoietin retards dna breakdown and prevents programmed death in erythroid progenitor cells}, Science \textbf{248} (1990), no.~4953, 378--381. \bibitem{Lajtha1959} L. G. Lajtha; \emph{On dna labeling in the study of the dynamics of bone marrow cell populations}, The Kinetics of Cellular Proliferation, Grune and Stratton, New York (1959), 173--182. \bibitem{Mackey1978} Michael C. Mackey; \emph{Unified hypothesis for origin of aplastic-anemia and periodic hematopoiesis}, Blood \textbf{51} (1978), no.~5, 941--956. \bibitem{MackeyHaurieBelairBOOK2003} Michael C. Mackey, Caroline Haurie, Jacques B{\'e}lair; \emph{Cell replication and control}, Nonlinear dynamics in physiology and medicine, Interdiscip. Appl. Math., vol.~25, Springer, New York, 2003, pp.~233--269. \bibitem{MackeyRudnickiJMB1994} Michael C. Mackey, Ryszard Rudnicki; \emph{Global stability in a delayed partial differential equation describing cellular replication}, J. Math. Biol. \textbf{33} (1994), no.~1, 89--109. \bibitem{MackeyRudnickiJMB1999} Michael C. Mackey, Ryszard Rudnicki, \emph{A new criterion for the global stability of simultaneous cell replication and maturation processes}, J. Math. Biol. \textbf{38} (1999), no.~3, 195--219. \bibitem{MahaffyBelairMackeyJTB1998} Joseph~M. Mahaffy, Jacques B{\'e}lair, Michael~C. Mackey; \emph{Hematopoietic model with moving boundary condition and state dependent delay: Applications in erythropoiesis}, Journal of Theoretical Biology \textbf{190} (1998), no.~2, 135--146. \bibitem{PujoBernardMackeySJADS2005} Laurent Pujo-Menjouet, Samuel Bernard, and Michael~C. Mackey, \emph{Long period oscillations in a {$G_0$} model of hematopoietic stem cells}, SIAM J. Appl. Dyn. Syst. \textbf{4} (2005), no.~2, 312--332. \MR{2173530 (2006e:34163)} \bibitem{PujoMackeyCRB2004} Laurent Pujo-Menjouet, Michael~C. Mackey; \emph{Contribution to the study of periodic chronic myelogenous leukemia}, Comptes Rendus Biologies \textbf{327} (2004), no.~3, 235--244. \bibitem{ReyMackeyCAMQ1993} Alejandro D. Rey, Michael C. Mackey; \emph{Multistability and boundary layer development in a transport equation with delayed arguments.}, Can. Appl. Math. Q. \textbf{1} (1993), no.~1, 61--81 (English). \bibitem{TanimukaiKimSakAlEH1997} S. Tanimukai, T. Kimura, H. Sakabe, et al.; \emph{Recombinant human c-mpl ligand (thrombopoietin) not only acts on megakaryocyte progenitors, but also on erythroid and multipotential progenitors in vitro}, Experimental Hematology \textbf{25} (1997), 1025--1033. \end{thebibliography} \end{document}