Stability analysis of a state-dependent delay differential equation for cell maturation: analytical and numerical methods.

We consider a mathematical model describing the maturation process of stem cells up to fully mature cells. The model is formulated as a differential equation with state-dependent delay, where maturity is described as a continuous variable. The maturation rate of cells may be regulated by the amount of mature cells and, moreover, it may depend on cell maturity: we investigate how the stability of equilibria is affected by the choice of the maturation rate. We show that the principle of linearised stability holds for this model, and develop some analytical methods for the investigation of characteristic equations for fixed delays. For a general maturation rate we resort to numerical methods and we extend the pseudospectral discretisation technique to approximate the state-dependent delay equation with a system of ordinary differential equations. This is the first application of the technique to nonlinear state-dependent delay equations, and currently the only method available for studying the stability of equilibria by means of established software packages for bifurcation analysis. The numerical method is validated on some cases when the maturation rate is independent of maturity and the model can be reformulated as a fixed-delay equation via a suitable time transformation. We exploit the analytical and numerical methods to investigate the stability boundary in parameter planes. Our study shows some drastic qualitative changes in the stability boundary under assumptions on the model parameters, which may have important biological implications.

Introduction

Maturation of cells, from the stem cell phase up to the fully mature phase, is an essential, and yet not completely understood, phenomenon. It is particularly interesting to understand what are the processes that promote a healthy equilibrium population state, or homeostasis, and what are those that may induce destabilisation of the equilibrium and appearance of periodic behaviour. We remark that oscillations play an important role in stem cell population dynamics in relation to hematological disorders like cyclic neutropenia and periodic myelogenous leukemia, see for instance Mackey (1978), Bernard et al. (2003), Pujo-Menjouet et al. (2005), and Adimy and Crauste (2012).

In this paper we investigate a stem cell model where maturation is described as a continuous process, and we show that the stability properties of the positive equilibrium are crucially affected by the choice of the maturation speed. In particular, we prove that different choices of the maturation speed imply qualitative and quantitative changes in the regions of stability of the equilibrium in parameter planes. We achieve this by developing analytical and numerical methods for the analysis of equilibria and stability of structured population models.

More precisely, the model is motivated by the blood cell production system and in particular by the production of white blood cells. It describes the process of division of cells by mitosis, their self-renewal, the differentiation of stem cells into progenitor cells (also called simply progenitors), and a continuous maturation process of the progenitors up to the transition into the mature cell compartment, see also Fig. 1. The mature cell population regulates the processes of division, self-renewal, differentiation and maturation. The development of progenitors, and in particular their self-renewal and maturation, depends on their maturity, which is described by a deterministically developing one-dimensional variable. In this way the progenitor population has continuous maturity structure. The latter is useful when it is not possible to divide the cell population into a finite number of discrete subpopulations. The continuous variable may represent a measurable quantity that evolves during maturation. This can be for instance age, size, or, as mentioned by Potten and Loeffler (1990) (to which we also refer for a short review of the biological terms), position in the tissue or “the weight of a specific protein per cell”. An advantage of keeping the model general, without specifying the maturity variable, is that it potentially applies to different cell types and laboratory measurements.

Fig. 1

Schematic representation of the model: at time t, w(t) and v(t) denote the total amount of stem cells and mature cells, respectively; u(t, x) is the amount of progenitor cells with maturity . The processes are indicated in the figure

We note that a similar concept of continuous maturity structure in the context of blood cells is considered by Adimy et al. (2005). Hematopoietic stem cells are divided into a resting and a proliferating phase, and the duration of the proliferating phase is assumed to depend on the maturity at commitment. Adimy and Crauste (2009) introduced regulation of the processes by growth factors, which are for instance proteins like G-CSF or EPO, that are controlled by the amount of mature cells. Adimy et al. (2010) considered an age-structured model with fixed duration of the proliferating phase. We refer to Pujo-Menjouet (2016) for a recent review of models for blood cells. All the above-mentioned models focus on the stem cell phase, which they structure by maturity. Moreover, they assume that cells divide only at the end of the proliferating phase, and division is assumed to be symmetric: upon division, two identical cells are produced and enter immediately the resting phase.

We consider a model for cell maturation where stem cells are assumed to be unstructured, whereas the maturity structure is introduced in the progenitor compartment. Moreover, we include the case of asymmetric cell division: two daughter cells are not necessarily of the same type, but they may belong to different compartments (stem or progenitor). In particular, this is achieved with the concept of fraction of self-renewal of stem cells (see also s(v) in Table 1), representing the fraction of daughter cells that are themselves in the stem cell compartment. The complementary fraction of daughter cells is assumed to enter the progenitor compartment with the initial maturity level.

Table 1

Specifications for stem cell rates, from Marciniak-Czochra et al. (2009) and Getto and Marciniak-Czochra (2015)

Description	Function
Net growth rate	q(v)=2s(v)-1d_w(v)-μ_w
Rate of inflow into progenitors	γ(v)=21-s(v)d_w(v)
Division rate	d_w(v)=p1+k_pv
Fraction of self-renewal	s(v)=a1+k_av

The parameter is the stem cell mortality rate. The parameters are nonnegative, with

A model motivated by these assumptions, but including a finite number of compartments to describe the progenitor phase, was established by Marciniak-Czochra et al. (2009). The finite number of compartments can be replaced with continuous maturity structure, where a progenitor cell is assumed to either undergo maturation or divide and produce daughter cells at the same maturity stage. In the literature, this has been done in two ways (see also Sect. 2 for more details). The first approach is by means of a partial differential equation (PDE) of transport type for the progenitor cells, coupled with two ordinary differential equations (ODE) for, respectively, stem and (fully) mature cells. This was first done by Doumic et al. (2011). Alternatively one can formulate the system as a delay differential equation (DDE) for the mature cells, in which the delay is state-dependent, coupled with the ordinary differential equation for the stem cells. A first formulation as a state-dependent delay differential equation (SD-DDE) was derived by Alarcón et al. (2011). Here we consider a slightly adapted SD-DDE formulation derived by Getto and Waurick (2016) (see again Sect. 2).

By guaranteeing differentiability of the functional inducing the DDE and an application of a theoretical result for SD-DDE of Walther (2003), Getto and Waurick (2016) showed that the model is well posed and that a linear variational equation can be associated with the solutions. The stem cell model has a trivial equilibrium and a positive equilibrium, which emerges from the trivial in a transcritical bifurcation. Recall that the principle of linearised stability allows to determine the stability of the equilibrium with respect to perturbed initial conditions by locating the roots of a characteristic equation with respect to the imaginary axis.

The cell maturation model considered in this paper is an instance of a physiologically structured population model in the sense of Metz and Diekmann (1986) and Diekmann et al. (2001). Such models are often formulated as Volterra functional equations coupled with DDE, in which the latter do not feature state-dependent delays. The principle of linearised stability in this formulation was proven by Diekmann et al. (2007) for the case of finite delay, and by Diekmann and Gyllenberg (2012) for the case of infinite delay. Diekmann and Gyllenberg (2007) gave conditions for differentiability of the right-hand side of a Volterra functional equation describing size-dependent cannibalism, thus providing the key to linearised stability analysis. A size-structured consumer–resource model was formulated by Diekmann et al. (2010, 2017) as a Volterra functional equation, and the linearisation computed analytically. For a special case of the model, the differentiability assumption required in Diekmann et al. (2007) was proven by Diekmann and Korvasová (2016). We stress that, apart from these papers, very few proofs of linearised stability exist for physiologically structured population models. Therefore, there are substantive grounds for developing methods for a more in-depth stability analysis. Here were develop both analytical and numerical methods.

Regarding analytical methods, in this paper we combine the linear variational equation derived by Getto and Waurick (2016) with theoretical results on linearised stability for SD-DDE by Hartung et al. (2006) and Stumpf (2016), to guarantee that the principle of linearised stability holds for the stem cell model in the DDE formulation. Then, we use the linear variational equation to derive characteristic equations for arbitrary, trivial and positive equilibria, respectively. We use these to derive the following results. The trivial equilibrium is stable in absence of the positive one and unstable in its presence, and the positive equilibrium is stable upon emergence through a transcritical bifurcation. This result is in agreement with the analysis of the multi-compartment model of Stiehl and Marciniak-Czochra (2011) and with the analysis of the PDE model of Doumic et al. (2011). Regarding the positive equilibrium, we also prove that in the right half-plane there are no roots of the characteristic equation outside of a compact set. These results can be combined with an argument of Diekmann et al. (1995) to conclude that roots can enter the right half-plane only through a compact subset of the imaginary axis, and thus the positive equilibrium can destabilise only if the latter occurs.

By the previous arguments, it makes sense to analyse the characteristic equation for the positive equilibrium on the imaginary axis, but the complexity of the characteristic equation in full generality motivates us to focus on special scenarios. Here, we consider a simplified characteristic equation where the dependence on maturity and on mature cell population of the progenitor development, given by maturation and self-renewal, can be neglected. In the population dynamical formulation, this simplification would transform the SD-DDE to a DDE with fixed delay, meaning a fixed positive time delay between leaving the stem cell compartment and entering the mature cell compartment. For this case, we show that the positive equilibrium can destabilise by a pair of conjugate roots crossing into the right half-plane, which gives evidence for a Hopf bifurcation.

To visualise the results (both analytical and, later on, numerical), we consider specifications of the model ingredients taken or adapted from the available literature. We single out pairs of parameters and we distinguish, in the associated parameter plane, the regions in which the positive equilibrium is stable and the regions in which it is unstable. These regions are shown via a curve, representing the stability boundary, that is parametrised by the position of the roots on the imaginary axis. This approach was presented for a prototype characteristic equation by Diekmann et al. (1995) and elaborated for a somewhat more general characteristic equation by Alarcón et al. (2014), and later by Diekmann et al. (2016). The characteristic equation analysed here is yet slightly more involved and we show some qualitative features that are not shown in the former references: depending on a third parameter, the curve describing the stability boundary has a unique minimum or is monotonically increasing.

For the analysis of the stability boundaries in more general cases, we propose numerical methods. The perhaps most generally applicable method for stability analysis in continuous population dynamics is the pseudospectral discretisation approach. The approach was first applied for studying the stability of the zero solution of linear DDE with fixed delay by Breda et al. (2005, 2013), see also Breda et al. (2015b) for a review, and applied to the linearisation of a size-structured consumer–resource model by Breda et al. (2015a). The pseudospectral discretisation was then applied by Breda et al. (2016a) to nonlinear delay equations, including both differential and integro-differential equations with fixed delay, and later extended to equations with infinite delay by Gyllenberg et al. (2018). With this method, the nonlinear DDE is approximated with a system of ODE, whose properties can be studied numerically with the package matcont for matlab, a well-established tool for the numerical bifurcation analysis of ODE developed by Dhooge et al. (2003, 2008). Here, we extend the approach to the case of SD-DDE. This is done by exploiting the uniform bound to the state-dependent delay, which ensures that the states of the dynamical system associated with the DDE are defined on a fixed and bounded interval. The latter interval is then discretised with a finite number of points and the state of the system is projected in the finite-dimensional space of polynomials via the pseudospectral approximation method. To the best of our knowledge, this is the only technique that allows the user to study numerically the stability boundaries of DDE with history-dependent delay using available software packages.

Since this paper contains the first application of the pseudospectral approach to SD-DDE, here we put a strong emphasis on numerical validation, remarking that this benefits the credibility of both the methods, numerical and analytical, involved here. Validation can be carried out for the case of constant maturation speed mentioned above (for which the analytical methods can also be applied), and for a more general case of maturation speed independent of the current value of maturity, that we are going to describe below.

Suppose that the speed of maturation of progenitors is still allowed to depend on the mature cell population, but it is independent of the maturity of the progenitor cell. This is a mathematically motivated simplification which may apply to a biological scenario in which the maturation rate of a progenitor cells at any stage is regulated only by external growth factors, which we assume regulated only by the population of mature cells. In the formulation of the population model, this means that the delay is still state-dependent, but the state-dependence is somewhat more explicit. In this case, a suitable time transformation allows to rewrite the SD-DDE as an equation with fixed delay.

We mention that, for a scalar equation with threshold-type delay, an analogue time transformation was first considered by Smith (1992, 1993). More recently, the transformation has been successfully applied to study the dynamical and bifurcation properties of delay models describing complex structured ecological systems, for instance by McCauley et al. (2008), Nelson et al. (2013), Bjørnstad et al. (2016) and Nisbet (1997). Apart from these applications, we could not find in the literature a rigorous statement of the correspondence of solutions of the two systems and of the related stability properties. Here we define the transformation explicitly introducing a parametrised map, and we show rigorously the invariance of equilibria and stability properties.

We then remark that DDE with merely fixed delays can be studied numerically with the package dde-biftool for matlab, developed by Engelborghs et al. (2002). Like matcont for ODE, dde-biftool allows the user to perform for instance the one-parameter continuation of equilibria and periodic solutions, the computation of eigenvalues and multipliers, and the continuation of bifurcation points in two parameters. Hence, in the case of maturation speed independent of maturity, there are three numerical methods for computing stability boundaries at disposal: pseudospectral method (and matcont analysis) applied to the SD-DDE, pseudospectral method applied after the transformation to a fixed-delay DDE, and dde-biftool after the transformation. We stress moreover that, in the case of a fixed delay in the untransformed model, the stability boundaries are available analytically and can be compared with the numerically computed curves. We use all of the above methods to compute stability boundaries and find that, upon changing the method, differences in the visualised boundaries are not distinguishable by eye. These tests support the validity of the pseudospectral approximation applied to SD-DDE, and justify the application of the method for studying the stability boundaries of more general instances of the model for which, as already mentioned, no other software package is available.

The stability boundaries presented here are the first to be computed for this model. They are carried out for rate specifications taken from the literature from Marciniak-Czochra et al. (2009), Doumic et al. (2011) and Stiehl et al. (2014). The specifications include different regulation mechanisms of the stem cell processes, like division or self-renewal, by the mature cell population. As for the progenitor maturation rate, for which we could not find an established formulation, we considered a new formulation that can incorporate different types of dependence on both maturity and mature cells. Of potential biological interest are some rather drastic changes in the curves upon changing these regulation modes. A cell biological interpretation of these changes and laboratory validation are some of the future perspectives of this research.

The paper is structured as follows. Section 2 introduces the model, see Eqs. (2.5)–(2.8). In Sect. 3 we recall some theoretical results about existence and uniqueness of solutions and the principle of linearised stability for equilibria. In Sect. 4 we compute the equilibria of the general model and the corresponding characteristic equations; we give conditions for the existence of a positive equilibrium which destabilises the trivial, and we show that the positive equilibrium is stable upon emergence. In Sect. 5 we restrict to a particular instance of Eqs. (2.5)–(2.8) with fixed delay; we give analytical formulas for the stability boundary of the positive equilibrium and we describe how the boundary changes qualitatively when the maturation rate changes. In Sect. 6 we consider the case of a maturation rate independent of maturity: we introduce a time transformation to define an equation with fixed delay with the same equilibria and the same stability properties. In Sect. 7, we explain how to approximate SD-DDE of the type (2.5)–(2.8) with the pseudospectral discretisation technique, and we specify the approximating system of ODE. Finally, in Sect. 8 we consider different assumptions on the maturation rate and we investigate numerically the effect on the stability boundary of the positive equilibrium. The proofs of the mathematical results of Sects. 3–6 are collected in the “Appendix”. We have tried to structure and highlight the paper such that each reader, inclined either to biological, computational or analytical aspects, should be able to find the relevant topics, with no need of going deeply into details about the other disciplines. In particular, we tried to present the main results and Sect. 9, which also contains a discussion of the potential biological relevance of the results, in a rather self-contained way.

The model

We consider a model for cell maturation studied in PDE formulation by Doumic et al. (2011) and in DDE formulation by Getto and Waurick (2016). A schematic illustration of the model is given in Fig. 1.

Cells are divided into stem cells, progenitor cells and fully mature cells. Stem cells and fully mature cells are unstructured in the sense that only their total amount, w(t) and v(t), respectively, affect the population dynamics. Progenitor cells, on the other hand, are structured by a one dimensional quantity denoted by x and taking on values in the finite interval . We let u(t, x) denote the density of progenitor cells with maturity x at time t.

We assume that the maturation process, including stem cell self-renewal and division, is subject to regulation by the mature cells v and by those only. In the terminology of Diekmann, Metz and collaborators [Metz and Diekmann (1986), Diekmann et al. (2001, 2003, submitted)], v plays the role of the environmental condition.

Let q(v) denotes the net growth rate of the stem cell population. This rate incorporates division, self-renewal, differentiation, which means transition to the progenitor compartment, and apoptosis, or mortality, of cells. The dynamics of the stem cells is given by the ordinary differential equationNext, let and d(x, v) be the maturation rate and net production rate, respectively, of progenitor cells. The rate d captures both self-renewal and decay of progenitor cells. The density u(t, x) of progenitor cells satisfies the PDEAt any time the total outflow from the stem cell compartment equals the inflow into the progenitor cell compartment. Therefore the PDE (2.2) must be supplemented by the following boundary condition:Finally, progenitor cells are assumed to become fully mature upon reaching maturity and this leads to the following equation for the fully mature cells:where is the per capita death rate of fully mature cells.

The system (2.1)–(2.4) specifies a physiologically structured population model. Observe that if v(t) is assumed to be a known function of time, the system (2.1)–(2.3) is a time-dependent linear system for w and u. The feedback described by (2.4) makes the full system into a nonlinear autonomous system.

The progenitor cell compartment can be eliminated from the system (2.1)–(2.4) by careful book-keeping. Notice that the progenitor cells reaching full maturity at at time t are the stem cells that differentiated and became progenitor cells with maturity at some time earlier plus those who have been born on the way due to self-renewal minus those who have died. Because the maturation rate g depends on the current density of fully mature cells, the maturation delay depends on the history of the fully mature cells, where we have used the standard notationif a function x is defined in .

To define we first define the function as the unique solution of the initial value problemgiven the history of the fully mature cells. Note that is the maturity of a progenitor cell s time units before it reaches full maturity at time t, that is, its maturity at time .

Because the maturation rate g is assumed to be positive, the function is monotone and therefore for any given history , the equationhas a unique solution , which is the time it took for a progenitor cell to mature from to , given that it reached at time t.

The book-keeping alluded to above amounts to integrating the PDE (2.2) along characteristics (Metz and Diekmann 1986). This was done by Getto and Waurick (2016). Substituting the resulting expression for into (2.4) they obtainedHere and in the rest of the paper, the notation denotes the derivative of a function with respect to its j-th argument.

The Eqs. (2.7) and (2.8) with y and given by (2.5) and (2.6), respectively, constitute the model which will be studied in the current paper.

Equations (2.7) and (2.8) consist of an ODE coupled with a differential equation with state-dependent and distributed delay (DDE) in which the state-dependent delay is defined via a threshold condition. We remark that a more explicit derivation of the DDE directly from first principles was attempted by the first author in a paper by Alarcón et al. (2011). However, that derivation erroneously missed the integral of in the exponent (the dilation factor, see Metz and Diekmann (1986)), as well as the ratio , which comes from the equation of fluxes at the boundaries of the maturation interval, taking into account that the maturation speeds differ at the boundaries and .

For a more in-depth analysis we will use the specifications for q and as given in Table 1. These ingredients were derived by Marciniak-Czochra et al. (2009) for a multi-compartment model describing hematopoietic stem cells producing leukocytes, and later considered by Getto and Marciniak-Czochra (2015) for the multi-compartment model as well as for the present model. They are based on the assumption that the individual stem cell division rate () and the fraction of self-renewal (s(v)) are regulated by a single external feedback mechanism through some signalling molecules called “cytokines”. The signal may either enhance the division rate or the fraction of self-renewal. The mechanistic derivation of the signal intensity (which justifies the form of the rates s(v) and in Table 1) was obtained by Marciniak-Czochra et al. (2009) through a quasi-steady state approximation. In this paper, we either consider a generic progenitor production rate d or neglect progenitor production by assuming . For g we will consider various choices that will be given below.

In the model of Marciniak-Czochra et al. (2009) stem cells are still described by (2.7), whereas progenitors are divided into a discrete number of compartments in which cells can divide, self-renew and differentiate (similarly as stem cells). Accordingly, the rates of inflow and net production of progenitor cells are described by rates with the same algebraic construction as the rate of inflow and net growth rate in Table 1, but possibly for different parameter values. Alarcón et al. (2011) and Doumic et al. (2011) replaced the discrete number of progenitor compartments by one progenitor compartment featuring continuous maturity structure as in (2.5) and (2.6). The dynamics within this compartment can be described by a classical transport equation. If the rates are specified as in Table 1, the model assumptions are equivalent to those of Doumic et al. (2011), and the results can be compared to this paper.

Linearised stability theorems

We show that the model can be analysed using results for a general class of DDE of the formwhere is defined on a suitable open subset for an appropriately chosen . The model (2.5)–(2.8) can obviously be written as an equation of the above type if we set and definefor , that will be specified below. We will consider two equilibrium states and the long-term behaviour of solutions starting closeby. To this aim we present two results, which refer to the class of DDE (3.1) and to the stem cell model (2.5)–(2.8), respectively. These state that, under certain conditions, the well-known principle of linearised stability holds. The proofs, along with associated results on existence and uniqueness of solutions for the respective systems, follow from a compilation of results from the literature that will be given in the following. We refer to this literature also for the mathematical details that we do not present here. Following Hartung et al. (2006), we say that f satisfies property (S) if the following hold:Moreover, we additionally introduce for f the following property (whose name stands for “strongly Lipschitz, uniformly on bounded sets”)The latter property is used for instance by Walther (2003). We define a function by introducing the notation

f is continuously differentiable;

each derivative extends to a linear map defined on ;

the map is continuous.

for any bounded set there exists some such that

Theorem 3.1

Suppose that f satisfies conditions (S) and (sLb). Then the DDE (3.1) is well posed (in the sense of Hartung et al. (2006)). If additionally is a constant function satisfying , then is locally asymptotically stable if all the roots of the characteristic equationhave negative real part, and unstable if the equation has a root with positive real part.

The proof essentially follows from the works of Diekmann et al. (1995), Walther (2003), Hartung et al. (2006), Stumpf (2016), and Getto and Waurick (2016). In the “Appendix”, in “Proof of Theorem 3.1”, we present a more precise argumentation.

We now introduce some assumptions on the model ingredients that are used to prove the principle of linearised stability for the model (2.5)–(2.8). We denote open balls by for some and some . First, we assume that g satisfies property (G) by Getto and Waurick (2016), that we recall here for convenience: there exist and open intervals I, J, with , such thatWith these assumptions, h can be defined as .

and is ;

is bounded on ;

on and .

Finally, we make the following additional assumptions on the model ingredients:Under these conditions, Getto and Waurick (2016) showed (see Theorem 1.13) that F satisfies (S) and (sLb). Hence Theorem 3.1 can be applied, and this immediately yields the following result.

are , bounded on , and

, and , are bounded on , whenever is bounded;

and are continuously differentiable, and Lipschitz on bounded sets.

Corollary 3.2

Suppose that g satisfies property (G), is a positive parameter, and g, d, and q satisfy assumptions (a)–(c). Then the stem cell model (2.5)–(2.8) is well posed (in the sense of Hartung et al. (2006)) and the solutions exist for all times. If additionally (w, v) is a constant function satisfying , then (w, v) is locally asymptotically stable if all the roots ofhave negative real part, and unstable if the equation has a root with positive real part.

For completeness, we repeat here some results from Getto and Waurick (2016) that will be useful later on. For , the setis a -submanifold of . As clear from the definition, X contains the segments of any solution after a certain finite time, including the constant functions satisfying the equilibrium conditions. For this reason X is called solution manifold. From (S) and (sLb) it follows that the delay functional is well defined by (2.5) and (2.6) and satisfiesfor all . The system (2.5)–(2.8) supplemented with the initial conditionis well posed, and a differentiable semiflow can be associated. Moreover, Getto and Waurick (2016) computed an expression for the derivative DF, which we will use in the present manuscript to analyse (3.3).

Equilibrium and stability analysis

Equilibria

Let us denote by (w, v) an equilibrium solution. If there is exactly one equilibrium, the trivial equilibrium . If the equilibrium conditions areWe assume that, on , the rate q decreases monotonically to a negative value. Then there exists a unique positive equilibrium if and only if . Thus, if we interpret q(0) as a bifurcation parameter, there is a transcritical bifurcation at .

Characteristic equations

We now present some results that will be used to compute and analyse the characteristic equation. All the proofs are collected in the “Appendix”. Given a function f defined on , we denote by its Laplace transform, defined byLet us denote by (w, v) the positive equilibrium. We will also use the equilibrium conditions (4.1), the latter of which allows to eliminate w.

Lemma 4.1

For the zero equilibrium the characteristic equation readsIf , for the positive equilibrium the characteristic equation iswhere

Exchange of stability at the transcritical bifurcation and a priori bounds

We first conclude that the positive equilibrium destabilises the trivial one.

Corollary 4.2

The roots of the characteristic equation for the trivial equilibrium are q(0) and . Thus the trivial equilibrium is locally asymptotically stable if and unstable if .

We next show that, if the positive equilibrium is sufficiently small, then the characteristic equation has exactly one root and this root is negative.

Lemma 4.3

Suppose that and consider (4.2). Then , f is analytic, and there exists a -function , defined in a neighbourhood of zero, with , such that, in a neighbourhood of (0, 0), (v, z(v)) is the unique solution. Moreover

Hence, by our earlier assumptions on q, if q(0) is sufficiently small then v is small enough for the corresponding unique root z(v) to be negative. Hence, the positive equilibrium is stable upon emergence. In summary we have shown an exchange of stability at the transcritical bifurcation.

We have shown that for q(0) sufficiently small the characteristic equation has no roots in the right half-plane. Our next question is whether the positive equilibrium can destabilise if q(0) moves further away from zero. In the following result we prove some a priori bounds for the roots. The result will be used to establish a maximal region of stability in parameter spaces.

Lemma 4.4

Suppose that and consider (4.2). Then there exists a K such that, for every root z with , one has .

We conclude this section by presenting the specification of the characteristic equation in some particular cases. If , we omit the first argument in the notation of g and write g(v) and instead of g(x, v) and respectively. Further down we consider the case where g is a constant function, the value of which we also denote by g. We will use similar notation for . We then get the following result.

Lemma 4.5

(a) If , the characteristic equation becomeswith(b) If moreover , , it becomes(c) If moreover q and are as specified in Table 1, it becomes(d) If moreover for some parameter ( in Table 1), it becomes(e) If moreover is as specified in Table 1, it becomes

Note that, after the specifications, the existence condition for the positive equilibrium can be translated toThen the existence region in the -plane is the region above the line . Upon crossing the line from below, a transcritical bifurcation with exchange of stability occur, see Fig. 4.

Fig. 4

Existence and stability of equilibria in the plane , for the rate specifications of Lemma 4.5(e), with and . The positive equilibrium exists for (no positive equilibrium in the striped region). In the large panel, the solid curves are the analytical curves (5.2) and (5.3) and show how the stability boundary changes qualitatively with : the positive equilibrium is stable below the curves, unstable above. The three upper panels contain some zooms of the curves: the instability region is shaded and the black dots are the numerical approximations of the curves computed by numerical continuation with the software dde-biftool. The rates correspond to the specifications (s) and (pv) in Table 2

Determination and analysis of the stability boundary for the positive equilibrium for an example with fixed delay

Destabilisation of the positive equilibrium

In this section we restrict to the case with fixed delay specified by the conditions of Lemma 4.5(e), and we analyse more in detail the stability boundary of the positive equilibrium. Multiplying (4.9) by , the following is proven.

Lemma 5.1

If p and are free parameters and the remaining ones are fixed, (4.9) can be expressed aswhere we introduced the increasing functionsa scaled stem cell mortality , and a scaled complex variable .

For the characteristic Eq. (5.1) can be written as the two real equationsSubstituting and solving with respect to r and m, we getNote that we should take care of the singularities of (5.2) and (5.3): we do this in the following result, together with a more complete analysis of the functions defining the denominators. For a visualisation of the following result, see Figs. 2 and 3. In the (m, r)-plane displayed in Fig. 3, the existence boundary corresponds to the positive vertical axis and the existence region to the positive quadrant.

Fig. 2

Graphs of (solid) and (dashed) for different values of

Fig. 3

Graphs of for and . Existence boundary and region correspond to positive vertical axis and positive quadrant, respectively

Lemma 5.2

For and the following hold.

For and , one has For ,

and ; in every there exists a unique zero; ; and on ; in every there exists a unique zero; and and on ;

; on ; ; has a unique zero in each ; ; on , hence there are no zeros in these intervals. There exists a unique zero in every .

; has no zeros in , and exactly one zero in every , . Moreover, for every , one has , on , there exists a unique zero in , and on ;

for one has and has no zeros in , neither in , a unique zero in each and in each .

The function has exactly the singularities , where and, for , the are the ordered remaining zeros of either or as determined in (a). In particular, for , one has for and for .

For the function is increasing on .

It follows from Proposition XI 2.13 of Diekmann et al. (1995) that, when crossing the curve from right—in the sense of increasing —to left, two complex conjugate roots cross from the left to the right half-plane, provided that the determinant in the following lemma is negative.

Lemma 5.3

One has

We have shown that, in a right neighbourhood of the transcritical bifurcation point, the positive equilibrium is stable. Since in Lemma 4.4 we have established a priori bounds for the roots, we know that roots can enter the right half-plane only through a compact subset of the imaginary axis. Now note that the signs of denominator in (5.2) and numerator in (5.3) agree, and the signs of the denominator of (5.3) and the determinant of Lemma 5.3 are opposite. It follows that, as long as the curves in the (m, r)-plane are in the positive quadrant, the sign of the determinant is negative. Hence we can combine the previous arguments with considerations on the ordering of the curves, defined on the intervals bounded by the singularities, to conclude that the boundary between the maximal region of stability and the maximal region of instability is given by the curve , withSee Fig. 4 for a plot. We refer to Chapter XI of Diekmann et al. (1995) for further details on the argumentation above.

Table 2

Summary of the parameter sets used in the numerical tests

	Parameters	Description
(s)_s	k_a=1, k_p=0	Stem cell regulated self-renewal, unregulated division (see Table 1)
(s)_d	k_a=0, k_p=1	Stem cell unregulated self-renewal, regulated division (see Table 1)
(pv)_0	k_1=0, k_2=0	Progenitor maturation independent of v (see (8.2))
(pv)_1	k_1=1, k_2=0	Progenitor maturation decreasing in v (see (8.2))
(pv)_2	k_1=0, k_2=1	Progenitor maturation increasing in v (see (8.2))
(px)_1	p_u(x)=0.3	First type of progenitor x-dependence (see (8.2) and Table 3)
(px)_2	p_u(x)=p	Second type of progenitor x-dependence (see (8.2) and Table 3)

See also Table 1 for specifications of stem cell rates, and Eq. (8.2) and Table 3 for the progenitor maturation rate

Analysis of the stability boundary

In the following, we use and as a short-hand notation for the one-sided right and left limit of a function f, respectively, whenever they exist, i.e.,From Lemma 5.2 it becomes clear that, if ,and, if ,with , see Fig. 3. Note that, for , . The stability region is unbounded in m, whereas it is bounded in r if , and unbounded if . In the following we present a rather complete analysis of the gradient of the stability boundary. We refer to Fig. 4 for a visual summary.

Lemma 5.4

One has on for , and on for .

It follows that we can alternatively define the stability boundary via (5.2) and (5.3) as , with

Lemma 5.5

For the stability boundary first decreases, then increases. Here, the minimum is assumed in values if , if , and if . For the curve is increasing.

Note that for there is a turning point after the minimum. Moreover, for both cases and , the proofs of the previous results show (although Fig. 4 does not) that and that there is a turning point also before the minimum.

Transformation to fixed delay for g independent of maturity

In this section we restrict to the case of a maturation rate g independent of maturity (i.e., ), and we establish a relation between the solutions of the SD-DDE and the solutions of a special equation with fixed delay. When , Eq. (2.8) readsIn this case, the explicit expressionallows us to write Eq. (2.6) aswithNext to (2.7) and (6.1), we consider the differential equationsfor , where the maximum delay is fixed by the model parameters, see (6.4). The system (6.5) and (6.6) is provided with the initial conditionfor belonging to the solution manifold. We will explicitly construct a transformation of the independent variable that, given a solution of (2.7) and (6.1), provides a solution of (6.5) and (6.6), and vice versa.

Given an interval with and a continuously differentiable function f defined on J with range in the domain of g, we define the transformation such thatNote that and is with for .

Similarly, we define the transformation such thatNote that and is with for .

Lemma 6.1

Let v be defined in and u be defined in , both with range in the domain of g. Define and . Then and are invertible with and .

We are now ready to specify the relation between the solutions of the initial value problems.

Theorem 6.2

Let (w, v) be a solution of (2.7) and (6.1) with initial condition (3.6) defined on , and let , . Then defined byis a solution of (6.5)–(6.7) on .

Vice versa, given a solution of (6.5)–(6.7) on , define , and let defined on such thatThen, (w, v) defined on byand extended to such that , , is a solution of (2.7) and (6.1).

We stress that, while the solution of the initial value problem (6.5)–(6.7) is uniquely defined from (w, v), the converse is not true: the construction of a solution (w, v) defined on from a given solution in general involves an arbitrary extension on . Such extension is required in the mathematical analysis to define the state of the dynamical system on the interval . However, this should not cause any problem in applications, where usually the initial data for system (2.7) and (6.1) is fixed by observations or chosen arbitrarily for experiments. From a biological point of view, the transformation translates the temporal time-scale into the “physiological” time-scale of progenitor cells (i.e., we measure the maturity level of cells instead of their age).

We consider the transformed functional G, which is given asNote that if and only if . We now let (w, v) denote the positive equilibrium. We remark that results similar to the following hold for the trivial equilibrium. We will not go into details.

Proposition 6.3

For it holds that if and only if .

The result is a simple corollary of the following

Lemma 6.4

For and we have

Therefore, system (2.7) and (6.1) and system (6.5)–(6.6) have the same equilibria, with the same stability properties.

Numerical stability analysis: the pseudospectral approach for differential equations with state-dependent delay

In this section we adapt the pseudospectral discretisation approach, presented by Breda et al. (2016a) in the case of nonlinear delay equations with fixed delay, to the case when the delay depends on the state. The main idea of spectral methods is to project a given equation defined on an infinite-dimensional space into a subspace of finite dimension, see Gottlieb and Orszag (1977). Pseudospectral techniques are particular kinds of spectral methods where the approximating space is chosen as the space of polynomials of a fixed degree , and the projection is done through collocation: a finite number of conditions is imposed on a set of points, called collocation nodes.

Breda et al. (2016a) showed that, by means of pseudospectral techniques, it is possible to obtain a finite-dimensional system of ordinary differential equations that approximates the dynamical properties of the original delay equation with fixed delay. Here, we extend the method to equations where the delay is state-dependent by exploiting the uniform bound of the delay (3.5). Thanks to the latter, the pseudospectral approach can be applied on the fixed interval , and the convergence analysis that holds in the fixed-delay case remains true also in state-dependent case, as we will explain below. For the sake of presentation, we summarise the main steps of the method and write explicitly the approximating system of ODE for Eqs. (2.7) and (2.8). We refer to Breda et al. (2016a) for further details and for more general types of delay equations.

The existence and uniqueness of solutions of (2.7)–(2.8) with initial condition (3.6) implies that the semigroup of solution operators is well defined, strongly continuous, and its infinitesimal generator iswhere X is defined in (3.4), see Crandall and Liggett (1971). Note that the action of is linear, while its domain X is defined via a nonlinear condition. Then, (2.7)–(2.8) with initial condition (3.6) is equivalent to the abstract differential equationfor , with initial conditionThe problems are equivalent in the sense that, if (w, v) is a solution of (2.7)–(2.8) and (3.6), then is a solution of (7.1)–(7.2) with , for all , and vice versa.

Equation (7.1) is a differential equation in the space of continuous functions. In order to obtain a numerical approximation, we fix a positive integer , called discretisation index, and define as the space of -valued polynomials of degree M. Then we consider the finite-dimensional space . The projection of the abstract differential Eq. (7.1) in the space iswhere , and is called the approximating generator. The operator is defined by imposing some collocation conditions on a finite set of nodes, as explained for instance by Gottlieb and Orszag (1977) and Boyd (2001), in the following way. We introduce a mesh of collocation nodes, , such thatGiven values , their interpolating polynomial on is uniquely determined, and it can be represented in Lagrange form as a linear combination of the values aswhere the coefficientsform the basis of Lagrange polynomials, as explained for instance by Quarteroni et al. (2007) in Chapter 8. Note that , where is the Kronecker’s symbol. Note also that, by linearity, the derivative of at the nodes can be expressed aswhereare the elements of the differentiation matrix corresponding to the Lagrange polynomial basis. Note that the elements can be explicitly computed by established numerical routines, as explained by Trefethen (2000).

Following Breda et al. (2016a), for all we define as the pair of polynomials satisfyingwhere is the functional defining the right-hand side of the DDE, as defined at the beginning of Sect. 3. Note that is well defined because every M-degree polynomial is uniquely determined by independent conditions. By exploiting the representations (7.4) and (7.5) into (7.3), it is easy to obtain the following result.

Proposition 7.1

The pseudospectral approximation of system (2.7)–(2.8) readsfor , , whereThe polynomials and are the pseudospectral projections of w(t) and v(t), respectively.

The system (7.6) consists of ODE, and its properties can be studied with well-established software like matcont for matlab.

In practice, the interpolation nodes are chosen asIn the literature, these points are known as Chebyshev extremal nodes, and they are defined as the abscissas corresponding to the maxima and minima of Chebyshev orthogonal polynomials. Chebyshev nodes ensure that the interpolation process converges on functions that are at least Lipschitz continuous, as explained for instance by Davis (1975) and Quarteroni et al. (2007). Furthermore, the interpolation of infinitely differentiable functions converges with spectral accuracy, viz. the error decays faster than for any integer k, as shown by Trefethen (2000, 2013).

The one-to-one correspondence of the equilibria of (2.7)–(2.8) and (7.6) can be proven in the same way as done by Breda et al. (2016a). Moreover, the fact that the states of the dynamical system are defined in allows us to use convergence results about interpolation on a bounded and fixed interval. Since the eigenfunctions associated with the linearised system are exponential, we can exploit the spectral convergence of the interpolation process and use the argument of Breda et al. (2005) to ensure the convergence of the associated characteristic roots. In particular, the eigenvalues corresponding to the linearisation of (7.6) approximate the rightmost eigenvalues of the linearisation of (7.1) with spectral accuracy, as proven by Breda et al. (2005, 2015a, b). In conclusion, from Theorem 2.4 and Corollary 2.7 of Breda et al. (2016a), by using results on interpolation on , the following result follows.

Proposition 7.2

If (w, v) is an equilibrium of (2.7)–(2.8), thenwithis an equilibrium of (7.6). Vice versa, if (7.7) is an equilibrium of (7.6), then (7.8) holds and (w, v) is an equilibrium of (2.7–2.8). Moreover, the characteristic roots and the bifurcation points associated with the equilibria of (7.6) converge to those of the corresponding equilibria of (2.7)–(2.8) with spectral accuracy when M increases.

In the next section, we first validate the pseudospectral approximation technique on some instances of system (2.7)–(2.8) that can be transformed via the time transformation introduced in Sect. 6. In particular, the numerical output of the pseudospectral approximation (for and tolerance ) was compared with the output of dde-biftool on the fixed-delay reformulation. All the curves in Figs. 4, 5, 6 and 7 were computed both with dde-biftool and with the pseudospectral method, and all the respective curves are indistinguishable to the eye. Since dde-biftool is a well-established package for numerical bifurcation of delay equations, we take the agreement of the numerical output as an important evidence of the validity of the pseudospectral discretisation approach applied to a type of equations (viz., SD-DDE) to which it has never been applied before. After the validation of the method, we exploit the pseudospectral approach to approximate the stability boundary of system (2.7) and (2.8) in some cases when the time transformation is not applicable and there is no well-established software available for bifurcation analysis.

Fig. 5

Stability boundary for parameter set (pv) and (s), for different values of . The positive equilibrium exists for , and it is stable in the region below the stability boundary

Fig. 6

Stability boundary for parameter set (pv) and different values of . Left: stem cell regulation (s); right: (s). Stability region below the boundary

Fig. 7

Stability boundary for parameter set (pv) and different values of . Left: stem cell regulation (s), stability region to the right of the boundary; right: (s), stability region to the left of the boundary

Numerically computed stability boundaries

For the numerical results in this section we consider the stem cell rates specified in Table 1 and we fix the following parametersFor this choice of parameters, condition (4.10) for the existence of a positive equilibrium (w, v) reduces to . As proven in Sect. 4, the positive equilibrium is stable upon emergence. We consider two different regulation mechanisms of stem cell processes (cfr. also Tables 1 and 2): In the following we investigate numerically the stability region of the positive equilibrium by including different types of v- and x-dependence in the maturation rate of progenitor cells.

regulated self-renewal, unregulated division: , ,

unregulated self-renewal, regulated division: , .

Table 3

Specification of x-dependent rates for progenitor cells, see also Fig. 8

	p_u(x)	a_u(x)
(px)_1	p_u(x)=0.3	Constant:	a_u(x)=1/3
		Linear:	a_u(x)=2/3(1-x)
		Quadratic:	a_u(x)=0.1+2.8(x-0.5)^2
(px)_2	p_u(x)=p	Constant:	a_u(x)=0.9
		Linear:	a_u(x)=0.99-0.18x
		Quadratic:	a_u(x)=0.855-0.54(x-0.5)^2

Fig. 8

Plot of from Table 3: linear (dashed), quadratic (dotted) and constant mean value (solid). Left: (px); right: (px)

In the modelling literature we could find very few specifications of how the maturation speed depends on maturity. Exceptions are the maturation rate proposed by Doumic et al. (2011), which is increasing, and by Alarcón et al. (2011), which is decreasing. Here we choose a maturation rate that allows to incorporate both cases. Hence, we specify the maturation rate of the formBy specifying the coefficients and we can vary the dependence on x; for illustration, we will consider two types of dependence that are mathematically convenient, see Table 3. We note that some parameter sets for the division rate and the fraction of self-renewal of cells in progenitor compartments are proposed by Stiehl et al. (2014).

In the following we also assume that the net production of progenitor cells is zero, i.e., , corresponding to a biological situation when progenitor death and self-renewal are in balance.

Constant progenitor maturation rate

We first consider the case of constant maturation rate of progenitor cells by taking (case (pv) in Table 2), and constant and , so that . This means that the maturation process is not regulated by mature cells, and the maturation speed is the same at all maturity levels. In this case, (2.7–2.8) reduces to an equation with one fixed discrete delay , whose bifurcation properties can be numerically analysed with the package dde-biftool.

The case (s) corresponds to the specifications of Lemma 4.5(e). The stability boundary was studied analytically in Sect. 5. The upper panels in Fig. 4 show the analytical stability boundary (solid line) and the numerically computed stability boundary obtained with dde-biftool (black dots) for different values of (we recall that for the choice (8.1)).

In the case (s), we studied the stability boundary for different values of varying between 0.1 and 50. Some examples are plotted in Fig. 5. The stability boundary moves towards higher values of p when the delay decreases. In this case the numerical simulations did not show any minimum: this behaviour is different from the case (s), where a minimum appears for small values of , see Fig. 4.

Maturation rate dependent on the amount of mature cells

We assume now that the maturation rate of a progenitor cell is regulated by the amount v of mature cells, but does not depend on the maturity of the cell itself: . Under these assumptions, we can apply the time transformation introduced in Sect. 6, so that the equilibria of system (2.7) and (2.8) and their stability properties are the same as in the fixed-delay system (6.5)–(6.6). We investigate numerically the latter system with dde-biftool. In particular, we take (8.2) withand we consider two cases (see also Table 2):

maturation rate decreasing with v: , ,

maturation rate increasing with v: , .

In the case (pv), the numerically computed stability boundary does not show a minimum for the parameter values considered in the simulations. Some examples are plotted in Fig. 6. When increasing , the boundary moves up. Interestingly, under assumption (s) we observe a change in the concavity of the boundary when is large enough, see Fig. 6 (right). As evident from the figure, by changing from concave to convex, the stability region of the equilibrium becomes much larger and perhaps unbounded in p, promoting stability of the system at equilibrium.

The type of regulation (pv) is elaborated by Doumic et al. (2011). In this case, the numerical results reveal that the stability region of the positive equilibrium has a qualitatively different shape compared to the cases examined so far, see Fig. 7. In both cases (s) and (s), the plots suggest the existence of an interval of values of such that the positive equilibrium is stable for every value of p, but in the first case this happens for large values of , while in the second case this happens for small values of . This motivated us to investigate also the case of simultaneous stem cell regulation (). No destabilisation was detected when varying the and p in the interval (0, 50). Thus, the regulation of the stem cell processes seem to have a significant impact on the shape of the stability region when the maturation rate is regulated by amount of mature cells only.

Maturation rate dependent on maturity and on the amount of mature cells

Finally, we investigate the case when the maturation rate depends explicitly on maturity x. We consider (8.2) with and maturity-dependent rates specified in Table 3. As clear from the specifications, we fix constant and we compare two different forms of the function (linear decreasing and quadratic) with their mean constant value, see also Fig. 8. A linear increasing function did not produce qualitatively different results, so we did not include the analysis here. Moreover, an increasing is not supported by the parameter sets proposed by Stiehl et al. (2014). We stress again that the transformation introduced in Sect. 6 does not apply to this case, and there is no well-established software for the numerical bifurcation analysis of (2.7)–(2.8). Therefore we study the system numerically by applying the pseudospectral discretisation method described in Sect. 7 and by studying the approximating system (7.3) with the numerical package matcont. In the following simulations, the discretisation index M is chosen as an integer in the interval [10, 15] that guarantees the convergence of the computation.

With parameter set (px), we considered different types of v-dependence corresponding to the parameter sets (pv), (pv) and (pv) summarised in Table 2. The results are summarised in Fig. 10. For parameter set (px), for illustration we only considered the case (pv) of maturation rate increasing with v (the case (pv) did not exhibit any more interesting behaviour). The output is plotted in Fig. 9.

Fig. 10

Stability boundary in the plane for case (px) in Table 3, when is constant (solid), linear (dashed) and quadratic (dotted). Different rows correspond to different types of v-dependence, see also Table 2. Stability region is below the boundary

Fig. 9

Stability boundary for case (px) in Table 3, when is constant (solid), linear (dashed) and quadratic (dotted), and v-dependence (pc). Left: stem cell regulation (s), stability region to the right of the boundary; right: (s), stability region to the left of the boundary

Introducing the dependence on maturity did not produce any qualitative change in the stability boundary compared to Figs. 5, 6 and 7, but affected the location of the boundary, sometimes in a relevant way, see for instance Fig. 9 (left). It is interesting to notice that the quadratic x-dependence in Fig. 10, (pv) (right), has a different effect (stability boundary moves up) compared to the other figures, where the stability boundary moves down. Also, in Fig. 9 we observe that the boundary associated with the quadratic dependence is located between those associated with linear and constant rates, differently from the outcomes of Fig. 10. These differences in the stability boundaries lead us to believe that the effects of maturity dependence are of great complexity and strongly dependent on the specifications of model parameters.

Discussion and outlook

In the present paper we have considered a model for the regulated maturation process of stem cells. We have developed analytical and numerical methods and used these to show that a stable positive equilibrium destabilises the trivial in a transcritical bifurcation, and that, upon further variation of parameters, the positive equilibrium may destabilise. Note that we could introduce a quantity R that can be interpreted as the expected lifetime net offspring production of a stem cell in absence of regulation. In the special case analysed in Sect. 5 (with the rates specified in Lemma 4.5(e), i.e., fixed delay only), this quantity would be and, by considering R as bifurcation parameter, the transcritical bifurcation occurs at as in classical population ecology and epidemiology, see for instance Diekmann et al. (1990). We note that, in the case of regular cell production, the healthy state at homeostasis corresponds to the positive equilibrium. If the model describes cancer stem cells, the positive equilibrium may represent a pathological situation, while the zero equilibrium may correspond to the absence of cancer tissue.

One motivation for the present research is the analysis of a six-compartment model of Marciniak-Czochra et al. (2009). By means of simulations the authors show that, under various choices of regulation mode and parameters, total cell numbers converge to a positive equilibrium, which is associated with homeostasis of the cell population. As remarked by Marciniak-Czochra et al. (2009), it is clear that convergence to a positive equilibrium is not feasible without regulation. In the case of two compartments, Getto et al. (2013) showed analytically that there exists a unique positive equilibrium in a certain region of the parameter space, while there is no positive equilibrium in the remaining region. Via a Lyapunov function it is proven that the positive equilibrium is globally stable under various regulation modes. Hence, Marciniak-Czochra et al. (2009) and Getto et al. (2013) found similar stability results, but they did not find destabilisation of the positive equilibrium (although Marciniak-Czochra et al. (2009) presented an example where the positive equilibrium is unstable on the whole existence region). These results are consistent with our findings in the sense that convergence to a positive equilibrium is possible when the stem cells population net growth rate is regulated.

Nakata et al. (2012) considered one unstructured progenitor compartment besides stem and mature cells, resulting in a three-compartment model. They proved that only in presence of the intermediate compartment the positive equilibrium can destabilise upon variation of parameters. Consistently with the analysis of Nakata et al. (2012), we have shown that adding a compartment between stem and mature cells can—depending on parameters—either preserve stability or allow for destabilisation. Such intermediate compartment is considered unstructured by Nakata et al. (2012), whereas it has a continuous maturity structure in the present paper. Comparing the stability boundaries of Nakata et al. (2012) with the ones found here, however, many substantial differences in shapes can be noted. It could be interesting to analyse these differences more precisely, in particular in relation with experimental data. Stiehl and Marciniak-Czochra (2011) and Nakata et al. (2012) showed that models with more than two unstructured progenitor compartments have, next to a zero and a positive equilibrium, also semi-trivial equilibria where the total amount of cells is positive, but where stem cells are absent. Apart from complicating the mathematical analysis, this feature may be questionable in terms of biological interpretation. As already observed by Doumic et al. (2011), models with a continuously-structured progenitor compartment do not allow for this type of semi-trivial equilibria and may therefore be more realistic in some biological contexts.

We remark that Doumic et al. (2011) studied a model similar to the one considered here, formulated as a transport equation. The paper contains results on boundedness of solutions, linearised stability, and persistence. Through formal linearisation (substituting exponential trial solutions) and numerical simulations, the authors found evidence for stability and for the destabilisation of the positive equilibrium via a Hopf bifurcation. In our stability analysis, we added a degree of reliability by quoting theorems that guarantee well-posedness and local correspondence of (in-)stability between original and linearised system, together with theorems that ensure that the model satisfies the necessary preconditions of the former theorems. Regarding the model ingredients, Doumic et al. (2011) focused mainly on a particular case motivated by the multi-compartment model,which is a special case of (8.2) with . Moreover, the progenitor production rate is not regulated by mature cells. Our results of stability on emergence hold for the more general model that incorporates a generic maturation rate g(x, v) and a progenitor production rate d(x, v) which may depend on mature cells.

One of the important results of Doumic et al. (2011) is the possibility of destabilisation of the positive equilibrium, which they accomplished by finding a parameter interval where destabilisation is possible, and by means of numerical examples. One of our aims here was to develop analytical and numerical methods for a quantification of the analysis, in particular the exact computation of destabilisation points and stability boundaries in parameter planes.

With analytical methods we have proven that the boundaries show interesting phenomena such as switches from boundedness to unboundedness in the stability region, and the appearance and disappearance of monotonicity upon variation of a third parameter, see Fig. 4. The presence of local minima corresponds to switches of the type stable-unstable-stable upon variation of a single parameter, as evident for instance from Fig. 4 (cases and ), when varying the parameter .

The numerical computations show a wider range of potentially interesting phenomena. For instance, the type of regulation of the maturation rate by mature cells has significant impact: note the drastic differences in the stability boundaries in Figs. 6 and 7. The pseudospectral discretisation method allowed us to conduct the numerical analysis for the fully general model, although, with our specifications of the rates, adding maturity dependence in the maturation rate seemed to affect the boundaries mainly quantitatively rather than qualitatively, see Figs. 9 and 10 .

The proven features regarding the shapes of the stability boundaries could be a good starting point to understand possible biological mechanisms relating parameters at the cell level with behaviours observable at the cell population level. The numerical methods would then allow for a wider biological analysis: it would be interesting especially to exploit the flexibility of the pseudospectral approximation to conduct further analyses with parameter estimates coming from the experimental data.

We remark that we could find very few references among the modelling literature for a mathematical description of the maturation speed. One possible reason for this may be the fact that numerical methods being able to handle maturation rates dependent on maturity are not widespread, and that computation times significantly increase for such rates. A second possible reason is the biological state of the art, since cell maturation is not completely understood yet. We hope that the exploration contained in this paper, showing that the maturation rate strongly affects the stability of the positive equilibrium, can contribute to raise interest in such biological questions.

On the mathematical side, it could be interesting to investigate if explicit representations of stability boundaries can be computed algebraically for the case of a non-constant maturation rate, and whether qualitative changes with respect to the constant case can be shown. The numerical computations, together with the fact that qualitative changes appear already when varying the value of the constant maturation rate, suggest an affirmative answer.

Moreover, there is a clear analytical and numerical evidence that the destabilisation of the positive equilibrium occurs through a Hopf bifurcation. It could be interesting to prove the Hopf bifurcation analytically. As a “warm-up”, for the fixed-delay case one could check whether the Hopf bifurcation theorem by Diekmann et al. (1995) is applicable. For the state-dependent case, one could study related Hopf bifurcation theorems by Hu and Wu (2010). We refer to Adimy et al. (2010) for a Hopf bifurcation analysis in the context of an SD-DDE modelling cell maturation, which uses a result from Eichmann (2006). An ongoing project is the proof of existence of periodic solutions as fixed points of the next-state operators for the stem cell SD-DDE.

For a progenitor maturation rate that depends on maturity and mature cells, at present, pseudospectral methods seem to be the only method of analysis. This paper represents the first application of the method to a SD-DDE. This was possible by exploiting the uniform bound of the delay (3.5), and by applying the discretisation technique used for fixed-delay equations to the interval . We argued that the convergence results proven in the fixed-delay case remain true in this new context. We also showed that, for fixed-delay equations, the stability curves obtained with pseudospectral methods are indistinguishable from the ones obtained with the well-established software package dde-biftoolEngelborghs et al. (2002). This agreement supports numerically the convergence of the pseudospectral approximations of the stability and bifurcation properties of equilibria.

After the successful investigation of the stability of equilibria, it would be interesting to exploit the features of the software packages for ODE to push the analysis forward and investigate also the periodic solutions emerging from the Hopf bifurcations, together with their stability and bifurcation properties. We recall that the proof of convergence of the approximation of periodic solutions is still under investigation, but the numerical tests of Breda et al. (2016a, b) support the conjecture of spectral accuracy for equations with fixed delay. We remark also that, when the ODE describing the evolution of the structuring variable cannot be solved explicitly, like in the maturity-dependent cases considered in this paper, the computation times may be significant. An interesting future perspective is the application of the pseudospectral approximation technique to complex structured models of the type proposed by Diekmann et al. (2010), with a special attention to the efficiency and computation times.

Finally, we remark that, in the spirit of physiologically structured population models, it could be interesting to reformulate the present model as a Volterra functional equation coupled with a DDE, for which the principle of linearised stability was proven by Diekmann et al. (2007) and Diekmann and Gyllenberg (2012). One of the advantages of that formulation would be a larger set of admissible initial conditions.