Transit and lifespan in neutrophil production: implications for drug intervention

TL;DR

This paper compares different mathematical models for neutrophil production, extending existing models with new parameters and techniques, and analyzes their stability and implications for drug response predictions.

Contribution

It introduces a generalized transit compartment model with decoupled maturation time and applies linear chain techniques to reformulate models as distributed delay DDEs.

Findings

Incorrect mean maturation time definitions in existing models

Reformulation as distributed delay DDEs enhances model analysis

Stability and bifurcation analyses inform pharmacological modeling

Abstract

We compare and contrast the transit compartment ordinary differential equation modelling approach with distributed and discrete delay differential equation models. We focus on Quartino's extension to the Friberg transit compartment model of myelosuppression, widely relied upon in the pharmaceutical sciences to predict the neutrophil response after chemotherapy, and on a QSP delay differential equation model of granulopoiesis. We extend the Quartino model by considering a general number of transit compartments and introduce an extra parameter which allows us to decouple the maturation time from the production rate of cells, and review the well established linear chain technique from the delay differential equation (DDE) literature which can be used to reformulate transit compartment models with constant transit rates as distributed delay DDEs. We perform a state-dependent time rescaling of the Quartino model in order to apply the linear chain technique and rewrite the Quartino model as a distributed delay DDE, which yields a discrete delay DDE model in a certain parameter limit. We then perform stability and bifurcation analyses on the models to situate such studies in a mathematical pharmacology context. We show that both the original Friberg and the Quartino extension model incorrectly define the mean maturation time, essentially treating the proliferative pool as an additional maturation compartment, which can have far reaching consequences on the development of future models of myelosuppression in PK/PD.

Tables

Table 1. Table 1 : Parameter values from Quartino2014 for the parameters of interest in this study. The mean value of the distributed delay in ( 2.31 ) is given by τ = n / a 𝜏 𝑛 𝑎 \tau=n/a , or τ = 106.41 𝜏 106.41 \tau=106.41 hours.

Parameter	Units	Typical estimate (% RSE)
$N_0$	109 cells/L	3.53 (5)
$k_{tr}=k_{\text{P}}$	h-1	0.03759
$\gamma$	-	0.444 (4)
$\beta$	-	0.234 (8)
$G_0$	ng/L	24.3 (8)
$k_e$	h-1	0.592 (32)
$k_{\text{ANC}}$	h/109 cells/L	5.64
$k_{\text{circ}}$	h-1	0.099
$k_{in}$	h-1	498.1792

Table 2. Table 2 : Hopf bifurcation points for varying γ 𝛾 \gamma , a 𝑎 a , and τ 𝜏 \tau for 𝐗 2 ∗ subscript superscript 𝐗 2 \mathbf{X}^{*}_{2} computed as in Figure  5 , Values highlighted in bold correspond to periods characteristic of patients with cyclic neutropenia.

$\gamma$ (-)	$a$ (hours-1)	$\tau$ (hours)	Period (days)
0.87851	0.03831	104.4	$18.00$
0.86766	0.03759	106.41	$18.32$
$0.78911$	$0.03243$	$123.33$	$21.0$

Table 3. Table 3 : Hopf bifurcation points for varying γ 𝛾 \gamma and n 𝑛 n for 𝐘 2 ∗ subscript superscript 𝐘 2 \mathbf{Y}^{*}_{2} computed as in Table  6 but considering the characteristic equation of the distributed DDE model ( 3.11 ). The other parameters were fixed at their homeostasis values given in Table  1 . The last row corresponds to the discrete delay case.

$\gamma$ (-)	$n$ (-)	Period (days)
4.61186	1	7.767
1.69292	1.5	13.53
1.21571	2	15.95
0.90252	3.5	18.08
$0.86766$	$\mathrm{𝟒}$	$18.32$
0.75312	10	18.98
0.72367	20	19.08
0.71461	30	19.09
0.71021	40	19.10
0.69754	“$\infty$”	19.11

Table 4. Table 4 : Hopf bifurcation points for varying γ 𝛾 \gamma and τ 𝜏 \tau for 𝐘 2 ∗ subscript superscript 𝐘 2 \mathbf{Y}^{*}_{2} computed as in Table  6 but considering the characteristic equation of the discrete DDE model ( 3.8 ). The other parameters were fixed at their homeostasis values given in Table  1 .

$\gamma$ (-)	$\tau$ (hours)	Period (days)
0.72631	99.8	18.00
$0.69754$	$106.41$	$19.11$
0.65534	117.75	21.00

Table 5. Table 5 : For each line the Hopf bifurcation point were computed changing the respective parameter and following the solution λ 𝜆 \lambda of Eq. ( 3.27 ) in the ( σ , ω ) 𝜎 𝜔 (\sigma,\omega) -plane numerically, the period was estimated by 2 ​ π / ω 2 𝜋 𝜔 2\pi/\omega and the parameters values at homeostasis were obtained from Craig2016c .

Parameter (units)	Homeostasis Value	Hopf Bifurcation	Period (days)
${\gamma }_Q$ (days-1)	$0.1$	$0.22791$	$36.69$
${\eta }_{N_P}^h$ (days-1)	$1.6647$	$7.4$	$60.31$
${\eta }_{N_P}^{min}$ (days-1)	$1.4060$	$0.81$	$25.86$
$V_{max}$ (-)	$7.8669$	$93$	$5.20$
$b_V$ (ng/mL)	$0.24610$	$0.0184$	$5.20$

Table 6. Table 6 : For each line the Hopf bifurcation point were computed changing the respective pair of parameters and following the solution λ 𝜆 \lambda of Eq. ( 3.27 ) in the ( σ , ω ) 𝜎 𝜔 (\sigma,\omega) -plane numerically, the period was estimated by 2 ​ π / ω 2 𝜋 𝜔 2\pi/\omega and the parameters values at homeostasis were obtained from Craig2016c .

Parameters  (units)	Hopf Bifurcation	Period (days)
$(V_{max},b_V)$  (-,ng/mL)	$(30,7.8)$	$5.20$
$(V_{max},b_V)$  (-,ng/mL)	$(61,16.1)$	$5.20$
$(b_{N_P},{\eta }_{N_P}^{min})$  (ng/mL,days-1)	$(0.065,1.1)$	$25.86$
$({\eta }_{N_P}^h,{\eta }_{N_P}^{\text{min}})$  (days-1,ng/mL)	$(2,1.1)$	$32.58$

Equations

$$\frac{\text{d}P}{\text{d}t}$$

$$\frac{\text{d}{T}_1}{\text{d}t}$$

$$\frac{\text{d}{T}_j}{\text{d}t}$$

$$\frac{\text{d}N}{\text{d}t}$$

$$\frac{\text{d}P}{\text{d}t}$$

$$\frac{\text{d}{T}_1}{\text{d}t}$$

$$\frac{\text{d}{T}_j}{\text{d}t}$$

$$\frac{\text{d}N}{\text{d}t}$$

$$\frac{\text{d}G}{\text{d}t}$$

$$\frac{\text{d}P}{\text{d}t}=\frac{\text{d}{T}_j}{\text{d}t}=\frac{\text{d}N}{\text{d}t}=\frac{\text{d}G}{\text{d}t}=0$$

$$k_{\text{P}}=k_{tr}.$$

$$k_{\text{P}}{P}_0=k_{tr}{P}_0=k_{\text{circ}}{N}_0.$$

$$\tau =\frac{n}{a}.$$

$$T_j=\frac{N_0{k}_{\text{circ}}}{a},j=1,\dots ,n,P_0=\frac{N_0{k}_{\text{circ}}}{k_{tr}},$$

$$k_{in}=G_0(k_e+k_{\text{ANC}}{N}_0).$$

$$\frac{\text{d}N}{\text{d}t}$$

$$=f\left(t,N(t),{\int }_0^{\infty }P(t-u)g_a^p(u)du\right).$$

$${\int }_0^{\infty }{g}_a^p(u)du=1.$$

$$g_a^p(t)=\frac{a^p{t}^{p-1}{e}^{-at}}{\Gamma (p)},$$

$$\tau =p/a,$$

$$\frac{\text{d}N}{\text{d}t}=f(t,N(t),P(t-\tau )).$$

$$\frac{\text{d}}{\text{d}t}{g}_a^p(t)$$

$$=a\left(\frac{a^{p-1}{t}^{p-2}{e}^{-at}}{\Gamma (p-1)}-\frac{a^p{t}^{p-1}{e}^{-at}}{\Gamma (p)}\right)$$

$$=a(g_a^{p-1}(t)-g_a^p(t)).$$

$$\frac{\text{d}}{\text{d}t}{g}_a^1(t)=\frac{\text{d}}{\text{d}t}(a{e}^{-at})=-a{g}_a^1(t).$$

$$T_j(t)={\int }_{-\infty }^{t}P(s)g_a^j(t-s)ds={\int }_0^{\infty }P(t-u)g_a^j(u)du,j=1,\dots ,n.$$

$$\frac{\text{d}N}{\text{d}t}=f\left(t,N(t),T_n(t)\right).$$

$$\frac{\text{d}{T}_j}{\text{d}t}$$

$$={\int }_{-\infty }^{t}P(s)a(g_a^{j-1}(t-s)-g_a^j(t-s))ds$$

$$=a(T_{j-1}(t)-T_j(t)),j=\{2,3,\dots ,n\}.$$

$$\begin{gathered} \frac{\text{d}{T}_1}{\text{d}t}=P(t)g_a^1(0)+{\int }_{-\infty }^{t}P(s)\frac{\text{d}}{\text{d}t}{g}_a^1(t-s)ds \\ =a(P(t)-T_1(t)). \end{gathered}$$

$$T_j(0)={\int }_{-\infty }^{0}P(s)g_a^j(-s)ds={\int }_0^{\infty }P(s)g_a^j(s)dsj=1,\dots ,n.$$

$$T_j(0)=P_0{\int }_0^{\infty }{g}_a^j(s)ds=P_0,j=1,\dots ,n.$$

$$\begin{array}[]{l}\displaystyle\frac{\textrm{d}P}{\textrm{d}t}=\left(k_{\textit{\tiny P}}(1-E_{\textit{Drug}})\left(\frac{N_{0}}{N(t)}\right)^{\!\gamma}-k_{tr}\right)P\\ \displaystyle\frac{\textrm{d}N}{\textrm{d}t}=-k_{\textit{circ}}N+k_{tr}\!\int_{-\infty}^{t}P(s)g_{a}^{n}(t-s)ds.\rule[-2.84526pt]{0.0pt}{22.76219pt}\end{array}$$

$$T_j(t)={\int }_{-\infty }^t\frac{k_{tr}}{a}P(s)g_a^j(t-s)ds,j=1,\dots ,n,$$

$$\frac{\text{d}N}{\text{d}t}=-k_{\text{circ}}N+a{\int }_{-\infty }^t\frac{k_{tr}}{a}P(s)g_a^n(t-s)ds,$$

Full text

We show that both the original Friberg and the Quartino extension model incorrectly define the mean maturation time, essentially treating the proliferative pool as an additional maturation compartment, which can have far reaching consequences on the development of future models of myelosuppression in PK/PD.

††journal: Journal of Pharmacokinetics and Pharmacodynamics

1 Introduction

In the pharmaceutical sciences, the concept of lag time, or the delay between the administration and the absorption of a drug, is a well-established phenomenon which is often accounted for Steimer1982 . Physiologically-based pharmacokinetic models incorporating absorption models like the ACAT or ADAM Agoram2001 ; Jamei2009 were indeed conceived and developed in part to account for the enterohepatic circulation that contributes to the delay in drug concentrations in the blood after oral administration. However, regardless of the administration of a xenobiotic, various forms of delays are present throughout physiological systems. In addition to pharmacokinetic lags, systems-level delays play an important role in determining the pharmacodynamic response to treatment. As examples, intracellular and intrinsic viral delays contribute to more complicated viral load decay in patients with human immunodeficiency virus being treated with antiretroviral drugs Dixit2004 , and the hematopoietic system displays multiple delays along the pathways from the pluripotent hematopoietic stem cells (HSCs) to terminally differentiated circulating cells Mackey1990 .

Granulopoiesis, the process of neutrophil production, in particular, exhibits multiple delays and has been studied in depth owing to the role neutrophils play in the innate (and adaptive) immune response Mantovani2011 . Within the pharmaceutical context, neutropenia is a toxic side effect of chemotherapy, and impacts heavily on treatment success and overall survival outcomes Sternberg2006 ; Gruber2011 . There is therefore an established interest in mathematical models that can predict the response to chemotherapeutic drugs Friberg2002 ; Foley2008 ; Schirm2014 and accurately represent the feedback mechanisms regulating neutrophil homeostasis Craig2016c ; Krinner2013 .

To maintain basal circulating neutrophil concentrations, multipotent progenitor HSCs in the bone marrow differentiate into the myeloid lineage on their way to becoming circulating neutrophils. After commitment, cells proliferate and undergo several divisions during a phase where cell numbers increase exponentially. After proliferation, neutrophil progenitors no longer divide. Instead, they grow in size and number of receptors before being sequestered into a marrow reservoir Rankin2010a , where they either die through apoptosis or transit into circulation Christopher2007 . Once they exit from the bone marrow, neutrophils circulate very transiently, with a half-removal time on the order of 7-10 hours VonVietinghoff2008 , as they either rapidly die or marginate into tissues Rankin2010a . Granulopoiesis is controlled by various cytokines, of which granulocyte colony-stimulating factor (G-CSF) is the principal actor Ward1999a . By binding to receptors on the neutrophil membranes, G-CSF regulates the rate at which neutrophils are released into circulation, and modulates up-stream factors (differentiation into the myeloid lineage, proliferation of upstream neutrophil progenitors, speed of maturation) to replenish and regulate the concentration of neutrophils in the bone marrow reservoir. G-CSF is then internalised by the neutrophils and removed from circulation. In the case of elevated circulating concentrations, G-CSF is also cleared via a linear, renal pathway Lyman2011 and these dual routes of elimination are important determinants of the PKs of G-CSF Craig2016c . An overview of the process of neutrophil production is given in Figure 1.

Mathematical representations of granulopoiesis (and other similar physiological delay systems) fall into three classes: transit compartment models where delays are represented via a chain of first-order ordinary differential equations (ODEs), distributed delay systems where integro-differential equations represent a delay that takes a range of values determined through some probability distribution Campbell2009 ; Hearn1998 ; Adimy2007 ; Vainstein2005 , or delay differential equation (DDE) systems where the present state depends on past states via fixed or state-dependent delays Foley2008 ; Brooks2012 ; Craig2016c (for more detailed discussions on the various models used in modelling hematopoiesis and chemotherapy-induced neutropenia, see Pujo-Menjouet2016 and Craig2017 , respectively).

Here we focus on two models of granulopoiesis in particular: the Quartino model Quartino2014 and the Quantitative systems pharmacology (QSP) model of Craig Craig2016c . The Quartino transit compartment ODE model accounts for the effects and PKs of endogenous G-CSF and is an extension of the widely-used Friberg model Friberg2002 ; Friberg2003 , while the QSP granulopoiesis model of Craig2016c is a state-dependant delay DDE model that incorporates the concentrations of unbound G-CSF and G-CSF bound to its neutrophil receptors. We will show that the Quartino model Quartino2014 can be reformulated as a distributed delay DDE, which becomes a discrete-delay DDE in a certain parameter limit. This reformulation of the Quartino model leads to some additional insight on parameter choices and will lead us to generalise this model.

Since the maintenance of homeostasis or the pathogenic shift towards disease-states depend on the longterm behaviour of a given system’s steady states, stability is an integral concept in physiology. In what follows, we will study the stability of these three major granulopoiesis model-types (transit compartment, distributed and discrete delay, and QSP) by demonstrating the relationships and equivalencies between all three formalisms and analysing the resulting distributed delay model to provide a better understanding of the role model selection plays within a treatment context. Accordingly, we will discuss how these stability results can impact the incorporation and delineation of the effects of interindividual variability. We will also provide a historical context for the origins of transit compartment models from distributed delay models and DDEs.

This paper is divided as follows. We begin in Section 2.1 with an extension to the common ODE transit compartment model before introducing the more general distributed and discrete delay models in Section 2.2. Therein, we discuss the linear chain technique (Section 2.2.1) used to recover a transit compartment model from distributed delay systems. We then briefly introduce our previously published QSP model including endogenous G-CSF negative feedback in Section 2.3. The stability of the transit compartment/distributed/discrete delay models and the QSP model is analysed in Section 3 before we undertake bifurcation analyses (Section 4), which we discuss within the pharmaceutical sciences context in Section 5. We conclude by discussing our results in Section 6. Many of the proofs are provided in the appendices at the end.

2 Modelling granulopoiesis: three different approaches to handling delays

2.1 Transit compartment model with endogenous G-CSF

The Friberg model Friberg2002 is perhaps the most well-known model of myelosuppression after chemotherapy in the pharmaceutical sciences Craig2017 . Five compartments are used to represent the HSCs and early progenitors, circulating neutrophils, and the transit between the proliferative and circulative states. A feedback mechanism on the rate of proliferation determines the extent of myelosuppression of the chemotherapeutic agent. The model has been shown to generically represent a variety of chemotherapeutic drugs Friberg2003 and has been widely adopted in PK/PD studies of anti-cancer drugs. We write a generalised version of this model as

[TABLE]

which reduces to the Friberg model if we set $k_{\textit{\tiny P}}=k_{tr}=a$ and $n=3$. Here, $P$ is the concentration of proliferating progenitors, $T_{j}$ is the $j^{\text{th}}$ post-mitotic transit compartment, and $N$ is the circulating neutrophil concentration (all in units of 109 cells/L), while $k_{\textit{\tiny P}}$ is the rate of proliferation in the progenitor cell pool, $k_{tr}$ and $a$ are the transit rates between the maturation compartments, and $k_{\textit{circ}}$ is the rate of neutrophil exit from circulation (all in units of h*-1*).

An extension to the Friberg model, which we will refer to as the Quartino model, is presented in Quartino2014 and models the myelosuppressive effects of chemotherapy on progenitor and circulating neutrophils, the endogenous G-CSF response, and the effect of the administration of a glucocorticoid to induce a rapid increase in G-CSF. A model schematic is given in Figure 2. For our purposes, we can discount the administration of the glucocorticoid prior to chemotherapy and ignore the corresponding model terms. We write a generalised version of the model as

[TABLE]

where $G$ is the circulating G-CSF concentration (ng/L), $k_{\textit{\tiny ANC}}$ is the neutrophil-dependent rate of G-CSF elimination (h*-1*), $k_{e}$ is the G-CSF nonspecific elimination rate (h*-1*), $(G/G_{0})^{\gamma}$ is the feedback on the proliferation rate from circulating G-CSF concentrations, and $(G/G_{0})^{\beta}$ reflects the G-CSF feedback on the maturation rate. In most of the current work we do not consider the chemotherapeutic agent and set $E_{\textit{Drug}}=0$, unless otherwise stated.

We let $P_{0}$, $N_{0}$ and $G_{0}$ denote the homeostasis values of $P$, $N$ and $G$ respectively, obtained by setting

[TABLE]

in (2.2). In both the Friberg and Quartino models, it is a modelling assumption that

[TABLE]

The condition (2.4) is required in the Friberg model (2.1) to ensure that $N=N_{0}$ at homeostasis, and in the Quartino model (2.2) to ensure that $G=G_{0}$ at homeostasis. If $G_{0}$ were not the homeostasis value of $G$, it would be hard to justify the $(G/G_{0})^{\beta}$ terms appearing throughout the model, and the model ought to take a different form. Consequently we enforce the condition (2.4) throughout, and always assume that $k_{tr}=k_{\textit{\tiny P}}$ as in Quartino2014 .

To see why we generalise the model by including a new parameter $a$, note that at homeostasis the rate of production of proliferating cells, the rate that cells leave proliferation to enter the first transit compartment, the rate they leave the last transit compartment to enter circulation and the rate that they leave circulation must all be equal. In both models this results in

[TABLE]

The production rate in (2.5) is completely independent of the maturation time of the cells; provided cells both enter and leave maturation at the rate given by (2.5), changing the maturation time $\tau$ would only change the total number of cells that are in maturation (which is $\tau k_{\textit{\tiny P}}P_{0}$), but will not change the production rate in (2.5). We will see below that at homeostasis the maturation time $\tau$ for both models is given by

[TABLE]

Fixing $a=k_{tr}$ leads to two related modelling problems. First if we regard $a=k_{tr}$ as known then, since $n$ is an integer, equation (2.6) only allows for certain discrete values of the delay $\tau$. On the other hand, if as is more usual we suppose that $\tau$ is known then choosing $n$ an integer and imposing that $a=k_{tr}$ in (2.6) uniquely determines the value of $k_{tr}$ in (2.6), which in turn determines the production rate in (2.5). But we already noted that the production rate at homeostasis $k_{\textit{\tiny P}}P_{0}$ and the maturation time $\tau$ are independent.

In Quartino Quartino2014 a mean maturation time is defined by $MMT=(n+1)/k_{tr}$. Presumably the authors counted $n$ transit compartments plus one proliferation compartment. By showing the equivalence of the generalised Quartino model (2.2) to a distributed delay DDE in Section 2.2 we will find both the mean and variance of the delay, and show that even if $a=k_{tr}$ the correct formula for the mean maturation time should be $MMT=n/k_{tr}$, corresponding to (2.6), and not the formula used in Quartino Quartino2014 .

In the following sections we will consider general values of the parameters $a$ and $n$, but will take the values of the remaining parameters from Quartino2014 ; these values are tabulated in Table 1. To satisfy the homeostasis conditions (2.3) we obtain

[TABLE]

and the parameter constraint

[TABLE]

At homeostasis the total number of cells in the $n$ maturation compartments is $N_{0}k_{\textit{circ}}n/a$. Dividing this by the production rate given by (2.5) gives the average maturation time $\tau=n/a$ as stated in (2.6).

Notice that if $a=k_{tr}$ as in Quartino2014 then at steady state we have $T_{j}=P$ for all the transit compartments. In Quartino2014 the model (2.2) is considered with initial conditions at time $t=0$ equal to the steady-state values (which is natural for a chemotherapy study before the chemotherapeutic agent is administered), but we will consider the behaviour of the model for general non-negative initial conditions. Proof of the positivity of solutions to the Quartino system (2.2) can be found in Appendix B.1.

2.2 Gamma-distributed delay representation of the transit compartment granulopoiesis models

Distributed delay DDEs come in many varieties, but a reasonably general form is

[TABLE]

In simpler examples $P(t)\equiv N(t)$, but $P(t)$ can also be a separate variable defined by its own differential equation (as is the case in the granulopoiesis models considered in this work). The function $g_{a}^{p}(u)$ is a probability density with

[TABLE]

So, rather than the dynamics of $N(t)$ being determined by the current value of $P(t)$, the integral distributes the effect of $P$ across its previous values. In this work we will restrict attention to the gamma distribution, though other distributions do arise, in particular the uniform distribution. We write the probability density function $g_{a}^{p}$ of the gamma distribution as

[TABLE]

where $\Gamma(p)$ is the gamma function. When $n$ is a positive integer $\Gamma(n)=(n-1)!$, and the gamma function generalises the factorial function to real numbers $p$ with $\Gamma(p)=(p-1)\Gamma(p-1)$ for any $p>0$. The real positive parameters $a$ and $p$ determine the shape and rate of the distribution with the mean delay $\tau$ given by

[TABLE]

and standard deviation $\sigma^{2}=p/a^{2}$. If $p$ and $a$ are taken to infinity with their ratio $\tau$ held constant then the variance decreases to zero and the probability density function $g_{a}^{p}(t)$ becomes narrower and taller and approaches the $\delta$-function $\delta(t-\tau)$. In this limit the distributed delay DDE (2.9) reduces to a discrete delay DDE

[TABLE]

So discrete delay DDEs can be thought of as a limiting case of distributed delay DDEs. We will see below that when $p=n$ an integer, we can rewrite a gamma distributed DDE as an ODE, so gamma distributed DDEs provide a link between ODEs and discrete delay DDE models.

2.2.1 The Linear Chain Technique

The linear chain technique is used to convert some distributed delay differential equations (DDEs) into a corresponding system of ordinary differential equations (ODEs), or vice versa. The technique dates back at least to the work of Vogel in the 1960s Vogel63 ; Vogel65 , and first appears in the English literature in the work of MacDonald MacDonald78 ; MacDonald89 who called the method the linear chain trick. Most authors continue to use that name, but we prefer linear chain technique, because, as we will see, there is a true equivalency between the differential equation systems, and no trick is involved. It is usually more convenient to formulate problems as ODEs for numerical simulation, but sometimes more convenient to formulate them as DDEs for analysis. The linear chain technique is well-known and used in population biology and mathematical epidemiology, but is as yet not as well-known in other fields. The method has been independently rediscovered several times over the decades, being referred to as the fixed boxcartrain method by Goudriaan Goudriaan1986 , and recently used by Krzyzanski Krzyzanski2011 in a pharmaceutical sciences setting. There are several variants on this technique, and descriptions can be found in many places including Jacquez2002 ; MacDonald89 ; Smith2011 , but the simplest application is for a gamma distributed delay, for which we will detail the steps here.

The probability density function (2.11) has the property that for $p\neq 1$

[TABLE]

While for $p=1$

[TABLE]

Models of the form (2.9),(2.11) can in principle be considered for any real positive value of $p$, but in practice nearly all authors only consider $p=n$ a positive integer (one exception is Campbell2009 ), because then equations (2.14),(2.15) allow the distributed DDE to be reduced to an ODE. To do this let

[TABLE]

Then equation (2.9) can be rewritten as an ODE

[TABLE]

Differentiating (2.16), using Leibniz rule for $j>1$ (noting that $g_{a}^{j}(0)=0$ for $j>1$) we obtain

[TABLE]

While for $j=1$ (noting that $g_{a}^{1}(0)=a$)

[TABLE]

Together equations (2.17),(2.18),(2.19) redefine the (nonlinear) distributed delay DDE (2.9) as a system of $n+1$ ODEs. General DDEs can be posed as infinite dimensional dynamical systems, which introduces considerable mathematical difficulties, so being able to reduce some DDE models to finite-dimensional ODEs is mathematically very advantageous.

To complete the relationship between the distributed delay DDE (2.9) and the system of ODEs (2.17),(2.18),(2.19) we should take some care with the initial conditions. The distributed DDE (2.9) has infinite memory, and so to solve as an initial value problem from time $t=0$ we need to define a history function $P(t)$ for all $t\leqslant 0$, so that the right hand-side of (2.9) can be evaluated. With $P(t)$ so defined, for the DDE and ODE reduction to have equivalent solutions, by (2.16) the ODE must have initial conditions

[TABLE]

If it is assumed that $P(t)=P_{0}$, a constant for all $t\leqslant 0$ then, using (2.10), we see that (2.20) reduces to

[TABLE]

It is natural to ask if we can also go the other way; does a solution of the system of ODEs (2.17),(2.18),(2.19), define a solution of the distributed DDE (2.9)? It follows immediately from (2.21) that a solution of the ODE system with initial conditions $P(0)=T_{j}(0)$ for $j=1,\ldots,n$ does define a solution of (2.9). The equivalence for more general initial conditions for the ODE has also been established; in that case the ODE initial conditions define a finite number of constraints on the history function $P(t)$ for $t\leqslant 0$, which do not uniquely define $P(t)$, and the ODE defines a solution of the distributed DDE (2.9) for all choices of $P(t)$ that satisfy the constraints Cooke1982 ; MacDonald89 .

2.2.2 Gamma-distributed and discrete delay representations of transit compartment granulopoiesis models

The linear chain technique of Section 2.2.1 can be applied to establish the equivalence between transit compartment ODE models and corresponding distributed delay DDEs. Consider first the distributed DDE system

[TABLE]

We define $T_{j}(t)$ by

[TABLE]

which corresponds to (2.16) with $k_{tr}P(s)/a$ replacing $P(s)$. Writing the equation for $N(t)$ as

[TABLE]

and applying the linear chain technique of Section 2.2.1 we obtain the generalised Friberg transition compartment model of myelosuppression (2.1). Taking $a=k_{tr}$ and $n=3$ gives the Friberg model as stated in Friberg2002 , as has already been noted in Belair2015 .

While it is necessary to set $a=k_{tr}$ in (2.22) to recover the model as stated in Friberg2002 , equation (2.1) defines a transit compartment model for other values of $a$ also, and both the system of ODEs (2.1) and the distributed DDE (2.22) can be considered for general values $a>0$.

The extended Quartino endogenous G-CSF model Quartino2014 as stated in (2.2) and discussed in Section 2.1 cannot be stated simply as a distributed delay DDE via the linear chain technique. The maturation time in the Quartino model instead of being constant is state-dependent with the rate constants for the passage through each transit compartment given by

[TABLE]

which varies as $G(t)$ varies; it reduces to the same value as for the Friberg model only if $G(t)=G_{0}$. In contrast, the derivation of (2.14), which is essential in the linear chain technique, requires that the rate constant $a$ (and the power $p$) be constant, so to apply the linear chain technique the profile of the probability density function must remain constant and cannot vary with time or the solution. Thus, while it might be tempting to consider a distributed DDE of the form

[TABLE]

if we set $a=k_{tr}(G(t))/G_{0})^{\beta}$, then it is not possible to reduce this model to a system of ODEs because the derivation of (2.14) fails when $a$ is time-dependent. Instead we could consider the model

[TABLE]

and apply the linear chain technique with

[TABLE]

to obtain the transit compartment model

[TABLE]

which is similar to the Quartino model (2.2), but missing the $(G(t)/G_{0})^{\beta}$ factors in all the $T_{j}$ transit terms, and consequently does not model the effect of G-CSF on the maturation rate.

To write the Quartino model (2.2) as a distributed DDE, we first remove the state-dependency of the delays by rescaling time. Define a new time $\hat{t}(t)$ by

[TABLE]

By Theorem B.1 the right-hand side of (2.29) is strictly positive for $t>0$ so $\frac{\textrm{d}\hat{t}}{\textrm{d}t}>0$ and the new time variable $\hat{t}(t)$ is a strictly monotonic increasing function of $t$. Then we see that

[TABLE]

Strictly speaking we should define new variables $\widetilde{G}(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)=G(t)$, but following common practice we suppress the tildes and reuse the same variable names. Applying the same time-rescaling to all the equations we rewrite the Quartino model (2.2) as

[TABLE]

We refer to (2.30) as the time-rescaled Quartino model. Since the time rescaling satisfies $\hat{t}(0)=0$, the initial conditions for the Quartino model (2.2) at $t=0$ and the time-rescaled Quartino model (2.30) at time $\hat{t}=0$ are the same, and these two equations given equivalent solutions.

The time-rescaled Quartino model (2.30) has constant transition rates between the transit compartments, and consequently we can apply the linear chain technique to derive (2.30) from

[TABLE]

by letting

[TABLE]

To define an initial value problem for the distributed delay DDE (2.31) we need to specify $N(0)$, $G(0)$ and $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)$ for $\hat{t}\leqslant 0$. This in turn defines initial conditions for both the time-rescaled Quartino model (2.30) and the Quartino model (2.2) with $T_{j}(0)$ given by evaluating (2.32) with $\hat{t}=0$. If $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)$ is constant for $\hat{t}\leqslant 0$ then (2.32) implies that $T_{j}(0)=k_{tr}P(0)/a$, so there is an immediate equivalence between all three models for such initial conditions. Even if the Quartino model (2.2) were considered with different initial conditions, there is still a direct equivalence to the time-rescaled Quartino model (2.30), and as noted at the end of Section 2.2.1, also to the distributed DDE model (2.31). Consequently we have three equivalent forms of the same model, with a direct correspondence between the solutions of the differential equation systems (2.2) and (2.30) and (2.31).

Recalling (2.12) the mean value of the distributed delay in (2.31) is $\tau=n/a$. The time rescaling (2.29) is trivial at homeostasis when $G(t)=G_{0}$, so this also implies that the mean maturation delay is $\tau=n/a$ in the Quartino model (2.2) (and fact that we already derived by a different argument in (2.6)). Fixing $a=k_{tr}$ only allows a very granular control of the mean delay in the ODE model by varying the integer $n$. Mathematically it is more convenient to fix the delay $\tau>0$ and use $n$ and $a$ to control the shape of the distribution. For the distributed DDE model (2.31) we do not even need $n$ to be an integer. Recalling (2.10), in the limit as $n\to\infty$ and $a\to\infty$ with $\tau=n/a$ fixed, the distributed delay DDE (2.31) reduces to the discrete delay DDE

[TABLE]

We remark that in the discrete delay DDE (2.33) the delay $\tau$ is constant in the rescaled time-variable $\hat{t}$, just as the (same) mean delay $\tau=n/a$ is constant in the distributed delay DDE (2.31). In contrast the mean maturation time $\alpha(t)$ in the Quartino model (2.2) varies with $G(t)$ and satisfies

[TABLE]

where $\hat{t}(t)$ satisfies (2.29). If $G$ is held constant (but not necessarily equal to $G_{0}$), this gives a mean maturation time $\alpha$ in the Quartino model (2.2) of

[TABLE]

For the case of time-varying $G(t)$, the evolution of the mean maturation delay $\alpha(t)$ is defined by a differential equation (A.2), which we derive in Appendix A, where we also show the similarities between this state-dependency and the explicit state-dependency in the QSP model (2.34). But, in the current work, the time-rescaling equation (2.29) will be sufficient for our purposes.

Notice that while the derivation of (2.33) makes sense when considering the limit of the shape of the probability density functions $g_{n/\tau}^{n}(t)$ as they approach the $\delta$-function $\delta(t-\tau)$ when $n\to\infty$, it is problematical if one interprets the distributed DDE (2.31) via the ODE system (2.28) or (2.30) since then the limiting process would correspond to taking the number of compartments $n$ to infinity while increasing the rate constants $a$ to infinity also.

Since the system (2.30) corresponds to the Quartino model (2.2) with time rescaled by (2.26), positivity of solutions is guaranteed by Theorem B.1. The correspondence between the distributed delay system (2.31) and the system (2.30) ensures positivity of solutions of (2.31) for integer $n$ only, but actually positivity can be established for all real $n>0$. The proof of the positivity of solutions to (2.30) is given in Theorem B.2 in Appendix B.1.

2.3 A QSP model of granulopoiesis and its regulation by G-CSF

As previously mentioned, DDEs are frequently relied upon to model granulopoiesis given the delays inherent to hematopoiesis. Here we focus on the Quantitative systems pharmacology model of Craig2016c , which has been shown to account for the dynamics of neutrophil production and its negative feedback relationship with G-CSF–both bound to receptors on the surface of neutrophils and freely circulating–in a variety of scenarios. The model is written as

[TABLE]

Here $Q(t)$ is the concentration of HSCs ($10^{6}$ cells/kg), $N_{\!R}(t)$ the concentration of neutrophils in the bone marrow reservoir ($10^{9}$ cells/kg), $N(t)$ the concentration of circulating neutrophils ($10^{9}$ cells/kg), $G_{1}(t)$ the circulating G-CSF concentration (ng/mL), and $G_{2}(t)$ the bound G-CSF concentration (ng/mL). Here, and throughout, the superscript h denotes the homeostasis value of a quantity. The system (2.34) is subject to the initial conditions (ICs) and history functions

[TABLE]

where $\varphi_{1,2}(t)\in\mathcal{C}^{0}$ and

[TABLE]

$I_{G}(t)$ models the administration of exogenous G-CSF. As described in Craig2016c , the self-renewal $\beta(Q)$ and amplification factor $A_{Q}(t)$ of the HSCs are given by

[TABLE]

and the rate at which HSCs differentiate into neutrophil precursors is determined by the circulating concentration of G-CSF

[TABLE]

The rate at which the neutrophil progenitors proliferate is given by

[TABLE]

where $\tau_{N_{\!P}}$ days is the time it takes for proliferation. After exiting proliferation, cells mature with rate

[TABLE]

where the maximal age of maturing neutrophils is $a_{N_{\!M}}$. Given $V_{\!N_{\!M\!}}(G_{1}(t))$ depends on the circulating concentration of G-CSF, the time it takes neutrophils to mature satisfies

[TABLE]

and the total time for the process of granulopoiesis is then the sum of the time to completion of each process, given by

[TABLE]

Maturing neutrophils are assumed to be subject to a constant death rate $\gamma_{N_{\!M}}$, and their amplification factor $A_{N}(t)$ is given by the integral equation

[TABLE]

The fraction of G-CSF bound to neutrophil receptors given by

[TABLE]

regulates the rate with which cells exit the marrow reservoir as

[TABLE]

Mature neutrophils die from the marrow reservoir with rate $\gamma_{N_{\!R}}$. Cells that transit into circulation are removed with constant rate $\gamma_{N}$.

Proofs of the existence, uniqueness, positivity, and boundedness of solutions to (2.34) are provided in Appendices B.2 and C.

3 Stability Analysis

We now perform stability analyses of the different models derived in the last section to determine what parameter values, if any, will render an equilibrium point unstable, most frequently by having sustained oscillations about it. This is done using a well-established technique, namely linearising about this equilibrium and then calculating at which parameter values the ensuing characteristic equation will have roots with positive real parts. This same technique is traditionally applied to systems of ODEs in which case the characteristic equation is a polynomial.

In general, the characteristic equation associated with an arbitrary distribution is transcendental and possesses an infinite number of roots. As we shall see, the advantage of an integer-order gamma distribution is to yield a characteristic equation which is also a polynomial, reflecting the fact, mentioned in Section (2.2.1) that the gamma distribution yields a system of ODEs. In the context of comparing and contrasting the different models, we obtain a “continuity” result of sorts as we determine that approximation in distribution does lead to approximation in stability diagrams (see Belair2015 for a similar continuity argument).

3.1 Characteristic equations for the Quartino endogenous G-CSF Models

Consider first the generalised Quartino model (2.2). Let

[TABLE]

be the vector of solutions so that (2.2) can be rewritten in vector form as

[TABLE]

where $\mathbf{F}(\mathbf{X})$ represents the right hand side of the Quartino model (2.2). Let $\mathbf{X}^{*}$ be an equilibrium of the system (that is that $\mathbf{F}(\mathbf{X}^{*})=0$), then define $\mathbf{Z}=\mathbf{X}(t)-\mathbf{X}^{*}$ and let $\mathds{J}(\mathbf{X}^{*})$ be the Jacobian of (2.2) evaluated at $\mathbf{X}^{*}$ ($\mathds{J}(\mathbf{X})=d\mathbf{F}/d\mathbf{X}$). Then linearising about $\mathbf{X}^{*}$ we obtain

[TABLE]

where nonlinear terms of order $\mathcal{O}(\|\mathbf{Z}\|^{2})$ are neglected. Seeking a nontrivial exponential solution $\mathbf{Z}(t)=\mathbf{C}e^{\lambda t}$ of (3.3) with $\mathbf{C}\in\mathds{R}^{n+3}$, a vector of constants, and $\lambda\in\mathbb{C}$, we obtain the characteristic equation

[TABLE]

where $\mathds{I}\in\mathds{R}^{(n+3)\times(n+3)}$ is the identity matrix. Evaluating the determinant in (3.4) leads to the characteristic equation, which is stated as Eq. (D.2) in Appendix D.1. This gives a polynomial of degree $n+3$ in $\lambda$ for the Quartino model (2.2), and a polynomial of degree $7$ if we set $n=4$, as in Quartino2014 . A steady state is unstable if any of the roots of this polynomial have positive real part.

The characteristic polynomial for the time rescaled Quartino model (2.30) is also a polynomial in $\lambda$ of degree $n+3$, and actually has a simpler form than the characteristic polynomial for Quartino model (2.2). But to derive this characteristic polynomial it is convenient to first consider the characteristic functions of the discrete delay DDE (2.33) and the distributed DDE model (2.31).

Let $\mathbf{Y}(t):=(P(t),N(t),G(t))^{\intercal}$ denote the vector of solutions of the discrete delay DDE  (2.33), and $\mathbf{Y}_{\tau}:=\mathbf{Y}(t-\tau)$ be the vector of delayed solutions. Then we can rewrite the the DDE (2.33) as

[TABLE]

in vector form. Similar to the ODE case, let $\mathbf{F}(\mathbf{Y}^{*},\mathbf{Y}^{*})=\mathbf{0}$ be a generic steady state. Define the variables $\mathbf{Z}\coloneqq\mathbf{Y}-\mathbf{Y}^{*}$ and $\mathbf{Z}_{\tau}\coloneqq\mathbf{Y}_{\tau}-\mathbf{Y}^{*}$ and denote the linearisation matrices of (3.5) computed at $(\mathbf{Y},\mathbf{Y}_{\tau})=(\mathbf{Y}^{*},\mathbf{Y}^{*})$ by $\mathds{A}$ and $\mathds{B}$. Linearising (3.5) about $\mathbf{Y}^{*}$ and using the variables $\mathbf{Z}$ and $\mathbf{Z}_{\tau}$ yields

[TABLE]

The linearisation matrices $\mathds{A}$ and $\mathds{B}$ from (3.6) are calculated in Appendix D.3. Seeking a nontrivial exponential solution $\mathbf{Z}(t)=\mathbf{C}e^{\lambda t}$ for equation (3.6), with constant $\mathbf{C}\in\mathds{R}^{3}$ and $\lambda\in\mathbb{C}$, we obtain the characteristic equation

[TABLE]

where $\mathds{I}$ is the identity matrix. Evaluating the determinant in (3.7) gives the transcendental characteristic equation

[TABLE]

where the coefficients $a_{2}$, $a_{1}$, $a_{0}$ and $b$ are computed in Appendix D.3.

In general equation (3.8) has infinitely many roots, corresponding to the infinite dimensional nature of DDEs. The treatment of these equations is made tractable because although there can be infinitely many complex numbers $\lambda$ that satisfy (3.8), it is well known that for any real number $\sigma$ there can only be finitely many solutions $\lambda$ with $Re(\lambda)>\sigma$ (see for example Lemma 4.2 in Smith2011 ). To determine stability we need to ascertain whether all the roots have $Re(\lambda)<0$.

Comparing the discrete delay DDE (2.33) with the distributed delay DDE (2.31), we see that they differ in only one term. Thus the linearisation of the distributed delay DDE (2.31) follows exactly the steps taken for the discrete delay DDE (2.33). Then, following MacDonald MacDonald89 , the characteristic equation for the distributed DDE (2.31) corresponds to (3.8) with the term $e^{-\lambda\tau}$ replaced by

[TABLE]

where $G(\lambda)$ is the Laplace transform of the gamma probability density function, and hence we obtain

[TABLE]

where the coefficients $a_{j}$ and $b$ computed in Appendix D.3 are the same as those for (3.8). Notice that if $b=0$ then (3.8) and (3.10) both reduce to the same cubic polynomial.

If $n$ is an integer, equation (3.10) is the characteristic equation of both the distributed DDE (2.31) and the equivalent time-rescaled Quartino ODE model (2.30). In that case, for $b\neq 0$, equation (3.10) can be written as

[TABLE]

a polynomial of degree $n+3$, which can be used to determine the stability of the steady-states of these models. But since the time rescaling (2.29) is monotonic this will also determine the stability of the steady-states of the Quartino model (2.2). Characteristic equations which reduce to polynomials, such as (3.10) with $n$ an integer, are said to be reducible MacDonald89 .

The characteristic equation (3.10) is also valid for the distributed DDE (2.31) when $n>0$ is not an integer. For general irrational $n$, equation (3.10) can have infinitely many roots, as is the case for the discrete DDE (3.8). But, if $n=p/q$ is rational, where $p$ and $q$ are co-prime integers then we observe rather odd behaviour. For an integer $m$, suppose that $n$ is rational with $n=p/q\in(m,m+1)$, which implies that $p\in(qm,qm+q)$. Solutions of (3.10) with $b\neq 0$ then satisfy

[TABLE]

(though not all solutions of (3.12) will necessarily solve (3.10) if $q$ is even). Since (3.12) is a polynomial in $\lambda$ of degree $p+3q$, the discrete delay DDE (2.31) has at most $p+3q\in(q(m+3),q(m+4))$ characteristic values $\lambda$ which satisfy (3.10) when $n=p/q$ is rational.

It is natural to think of the discrete DDE (2.33) as the limit as $n\to\infty$ with $a=n/\tau$ of the distributed DDE (2.31), and indeed with $a=n/\tau$ we have

[TABLE]

so the characteristic equation (3.11) for the distributed DDE approaches the characteristic equation (3.8) of the discrete DDE as $n\to\infty$. However, if one considers $n$ varying across the real numbers this is not a smooth limit as $n$ transitions between the rationals and irrationals. Consequently, even though the distributed DDE model (2.31) is valid for general real $n$, most authors, even when considering the behaviour as $n$ is varied or as $n\to\infty$ mainly restrict attention to the case where $n$ is an integer Beretta16 ; Campbell2009 ; MacDonald89 .

3.2 Stability analysis for the Quartino endogenous G-CSF model and related forms

The generalised Quartino model (2.2), has two steady states. Assuming that $k_{tr}=k_{\textit{\tiny P}}$ as in (2.4) for the reasons already stated, and considering the model in the form (3.2) with vector solution $\mathbf{X}(t)=(P(t),T_{1}(t),\dots,T_{n}(t),N(t),G(t))\in\mathds{R}^{n}$ these are given by

[TABLE]

where $N_{0}$ is given by (2.8).

The time-rescaled Quartino model (2.30) has the same steady states $\mathbf{X}^{*}_{1}$ and $\mathbf{X}^{*}_{2}$, since a monotonic rescaling of time does not affect equilibria.

The discrete and distributed delay Quartino DDE models (2.31) and (2.33) have related equilibria, but in fewer space dimensions, since these models do not include transit compartments. In the $\mathbf{Y}(t)=(P(t),N(t),G(t))\in\mathds{R}^{3}$ notation of (3.5) these are given by

[TABLE]

If $n$ is a positive integer the distributed DDE (2.31) is equivalent to the Quartino model (2.2) and the steady states $\mathbf{Y}^{*}_{1}$ and $\mathbf{Y}^{*}_{2}$ correspond exactly to $\mathbf{X}^{*}_{1}$ and $\mathbf{X}^{*}_{2}$ as defined in (3.14) and (3.15) for the appropriate $n$, and with the values of $T_{j}$ following from (2.32). We have the following stability result for these equilibria.

Theorem 3.1

Provided the parameters satisfy the constraints (2.4) and (2.8)

For the distributed delay DDE (2.31) and the discrete delay DDE (2.33) the steady state $\mathbf{Y}^{*}_{1}$ is locally asymptotically stable if $\gamma<\beta$ and unstable if $\gamma>\beta$, and the steady state $\mathbf{Y}^{*}_{2}$ is unstable if $\gamma<\beta$. 2. 2.

For the Quartino model (2.2) and the time rescaled Quartino model (2.30) the steady state $\mathbf{X}^{*}_{1}$ is locally asymptotically stable if $\gamma<\beta$ and unstable if $\gamma>\beta$, and the steady state $\mathbf{X}^{*}_{2}$ is unstable if $\gamma<\beta$.

Proof

At $\mathbf{Y}^{*}_{1}$ from (D.10) we have $b=0$, thus from (3.8) and (3.10) the characteristic equation for both the discrete and distributed DDE models reduces to

[TABLE]

and the stability of $\mathbf{Y}^{*}_{1}$ is the same for both models. If $\gamma>\beta$ then from (D.10) we have $0>a_{0}=h_{1}(0)$, while $h_{1}(\lambda)\to+\infty$ as $\lambda\to+\infty$. Hence, by the intermediate value theorem there exists $\lambda>0$ such that $h_{1}(\lambda)=0$, and thus the steady state is unstable.

If $\gamma<\beta$ from (D.11) we have $a_{2}>0$, $a_{0}>0$ and $a_{2}a_{1}>a_{0}$ and it follows from the Routh-Hurwitz criteria MacDonald89 that $Re(\lambda)<0$ for all characteristic roots of (3.18), and hence $\mathbf{Y}^{*}_{1}$ is stable.

For the steady state $\mathbf{Y}^{*}_{2}$ when $\gamma<\beta$ equation (D.11) yields $a_{0}=0$, and $b>0$. For the discrete delay DDE (2.33) the characteristic equation reduces to

[TABLE]

Then $h_{2}^{\infty}(0)=-b<0$, while $h_{2}^{\infty}(\lambda)\to+\infty$ as $\lambda\to+\infty$, and the intermediate value theorem again implies that the steady state is unstable. For the distributed DDE (2.31) a similar proof shows instability using the characteristic function (3.11) becomes

[TABLE]

which has a positive leading coefficient and is negative when $\lambda=0$, so again the intermediate value theorem shows that the steady state is unstable.

The steady states $\mathbf{X}^{*}_{1}$ and $\mathbf{X}^{*}_{2}$ of the Quartino model (2.2) have the same stability as those of the time-rescaled Quartino model (2.30), as the monotonic time-rescaling does not change the stability, though it will change the values of the characteristic roots. But the time rescaled model (2.30) has its characteristic roots given by the degree $(n+3)$ polynomial (3.11) which has the same roots as the characteristic equation (3.10) of the distributed delay DDE, and hence $\mathbf{X}^{*}_{j}$ and $\mathbf{Y}^{*}_{j}$ have the same stability. ∎

For the standard parameters, as given in Table 1, we have $\gamma>\beta$ so Theorem 3.1 implies that the neutropenic steady states $\mathbf{X}^{*}_{1}$ and $\mathbf{Y}^{*}_{1}$ are unstable in all these models. Proving directly that the homeostatic steady states $\mathbf{X}^{*}_{2}$ and/or $\mathbf{Y}^{*}_{2}$ are stable when $\gamma>\beta$ is difficult, but we can compute the roots of the characteristic equations, and these are shown in Figure 3. We see that the homeostatic steady states are indeed stable when $\gamma>\beta$. Moreover the characteristic roots for the time-rescaled transit compartment Quartino model converge to the characteristic roots for the discrete delay DDE as $n$ increases. Although the steady state is stable for all the models, we see that it becomes less stable as $n$ increases, with the real part of the characteristic values tending to increase with $n$. The phenomenon of loss of stability for fixed delay $\tau$ as $n$ is increased has long been known, but remains an area of active interest MacDonald89 ; Campbell2009 ; Beretta16 .

In Section 4 we will study how increasing $n$ can make the system more susceptible to loss of stability through bifurcations. However, here we point out that the change in stability observed in Theorem 3.1 depending on whether $\gamma>\beta$ or $\gamma<\beta$ is not a bifurcation in the usual sense. When $\gamma=\beta$ the model is degenerate with the progenitor equation reducing to $\frac{\textrm{d}P}{\textrm{d}t}=0$.

In the proof of Theorem 3.1 we made use of the relationships between the different model formulations to greatly simplify the derivation of the stability results. In particular the simpler forms of the characteristic equation for the time-rescaled Quartino (2.30) and DDE models (2.31) and (2.33) makes these much easier to work with. It is nevertheless possible to directly derive stability results for the Quartino model (2.2) at least for the steady state $\mathbf{X}^{*}_{1}$, though the proofs are much more involved. We include those results in Appendix D.2 for completeness.

3.3 Stability of the QSP model of granulopoiesis

Here we perform the linear stability analysis of the steady states of the QSP model defined by the DDE system (2.34), without any exogenous G-CSF, i.e. $I_{G}(t)=0$. Similar to Section 3.2, let $\mathbf{X}(t)\coloneqq(Q(t),N_{\!R}(t),N(t),G_{1}(t),G_{2}(t))^{\intercal}$ be the vector solution of (2.34) and $\mathbf{X}_{\sigma}\coloneqq\mathbf{X}(t-\sigma)$ denote a vector of delayed solutions. Then the DDE system defining the QSP model (2.34) can be rewritten in vector form as

[TABLE]

Parameters changes to the model lead to different steady states in equation (3.19). Let the steady state computed at homeostasis be written as $\mathbf{X}^{h}\equiv(Q^{h},N_{\!R}^{h},N^{h},G_{1}^{h},G_{2}^{h})$, and denote a generic steady state by $\mathbf{X}^{*}\equiv(Q^{*},N_{\!R}^{*},N^{*},\!$ $G_{1}^{*},G_{2}^{*})$.

To linearize (3.19) around a steady state $\mathbf{X}^{*}$ we define the variables $\tau(t)=\tau_{N}(t)-t$ and $u=s+\tau(t)$ to rewrite the amplification factor (2.39) as

[TABLE]

Thus we approximate the amplification factor (3.20) through the linearisation of the proliferation function (2.37) given by

[TABLE]

where $\eta_{N_{\!P}}^{\prime}\equiv d\eta_{N_{\!P}}/dG_{1}$. Further, since it does not affect the local stability of the steady state Cooke1996 , we freeze the state-dependent delay at its steady state value

[TABLE]

Using (3.21) and (3.22) together with the distributed delay variable defined by

[TABLE]

equation (3.20) becomes

[TABLE]

As a consequence of the approximation in (3.23), we can rewrite (3.19) as

[TABLE]

where

[TABLE]

Let $\mathbf{X}^{*}$ be a generic steady state of (3.24), defined by $\mathbf{f}(\mathbf{X}^{*},\mathbf{X}^{*},\mathbf{X}^{*},\mathbf{X}^{*},\mathbf{X}^{*})\!$ $=\mathbf{0}$. Define the variables $\mathbf{Z}\coloneqq\mathbf{X}-\mathbf{X}^{*}$, $\mathbf{Z}_{\sigma}\coloneqq\mathbf{X}_{\sigma}-\mathbf{X}^{*}$ and $\tilde{\mathbf{Z}}\coloneqq\tilde{\mathbf{X}}-\mathbf{X}^{*}$ and denote the linearisation matrices of (3.24) with regards to $\mathbf{X}$, $\mathbf{X}_{\tau_{Q}}$, $\mathbf{X}_{\tau_{N}}$, $\mathbf{X}_{\tau_{N_{\!M\!}}}$, $\mathbf{X}_{\tau}$ and computed at $\mathbf{X}^{*}$, respectively by $\mathds{A}$, $\mathds{B}$, $\ldots$, $\mathds{E}$. Linearising (3.5) about $\mathbf{X}^{*}$ and using the variables $\mathbf{Z}$, $\mathbf{Z}_{\sigma}$ and $\tilde{\mathbf{Z}}$ yields

[TABLE]

The linearisation matrices $\mathds{A}$, $\mathds{B}$, $\ldots$, $\mathds{E}$ from (3.25) are computed in Appendix D.4. Seeking a nontrivial exponential solution $\mathbf{Z}(t)=\mathbf{C}e^{\lambda t}$ for equation (3.25), with constant $\mathbf{C}\in\mathds{R}^{5}$ and $\lambda\in\mathbb{C}$, we obtain the characteristic equation

[TABLE]

where $\mathds{I}$ is the identity matrix and

[TABLE]

For $\lambda=re^{i\theta}$, with $r,\theta\in\mathds{R}$, we have $\lim_{r\rightarrow 0}f(re^{i\theta})=1$. We rearrange equation (3.26) in the form $\det(\mathds{F})=0$. To calculate the matrix $\mathds{F}$, with terms $F_{ij}$ for $i,j=\{1,2,\ldots,5\}$, note that some terms of the linearisation matrices are symmetric while others are antisymmetric, namely $A_{43}=A_{42}$, $A_{52}=-A_{42}$, $A_{53}=-A_{42}$, $A_{35}=-A_{25}$, and $D_{24}=-A_{24}$ (see Appendix D.4). Using this fact, we can then compute the terms

[TABLE]

[TABLE]

and write the matrix $\mathds{F}$ as

[TABLE]

Defining the constants

[TABLE]

and the functions

[TABLE]

[TABLE]

the characteristic equation $\det(\mathds{F})=0$ becomes

[TABLE]

The solutions $\lambda\in\mathbb{C}$ to (3.27), the characteristic roots, determine the stability of the steady state $\mathbf{X}^{*}$ from (3.19)/(3.25). To evaluate the stability numerically, we write $\lambda=\sigma+i\omega$, with $\sigma\in\mathds{R}$ and $\omega\in\mathds{R}$, and then compute the roots of (3.27) in the $(\sigma,\omega)$-plane using the Matlab subroutine fsolve Mathworks . As illustrated in Figure 4, at homeostasis all the characteristic roots $\lambda$ have negative real part, and so the homeostatic steady state $\mathbf{X}^{h}=(Q^{h},N_{\!R}^{h},N^{h},G_{1}^{h},G_{2}^{h})$, defined by $\mathbf{F}(\mathbf{X}^{h},\mathbf{X}^{h},\mathbf{X}^{h},\mathbf{X}^{h})=\mathbf{0}$ in (3.19) is locally asymptotically stable.

4 Bifurcation studies

Bifurcation analysis, or the study of the qualitative changes to the behaviour of a system given a change to parameter values, is a fundamental dynamical systems concept Meiss2007 . Accordingly, studying bifurcation points can be a powerful tool in the life sciences to shed light on underlying parameter relationships and better understand the robustness of a system with regards to stability.

Historically, bifurcation analysis has been applied to study hematological pathologies and has provided valuable insight into the origins of disorders like cyclic neutropenia, a disease associated with dangerously low neutrophil counts and mouth blistering Dale2015 where a patient’s ANCs oscillate with a period of around 21 days. These oscillations have been shown to correspond to a periodic orbit that appears through a loss of stability after the system undergoes a Hopf bifurcation Colijn:2 ; Foley2009b . In the following sections, we perform bifurcation analyses on the equivalent forms of the Quartino endogenous G-CSF model and the QSP granulopoiesis model (2.34) to ascertain how changing parameter values modify the stability of each system, giving insight into the potential effects of PK variability on a physiological system and helping to understand pathophysiology of diseases.

4.1 Bifurcation in the equivalent expressions of the Quartino endogenous G-CSF model

We begin by investigating whether parameter changes in the equivalent expressions of the Quartino endogenous G-CSF model can lead their steady states $\mathbf{X}^{*}_{2}$ and $\mathbf{Y}^{*}_{2}$ given by (3.15) and (3.17), respectively, to lose stability. For this, we let $\lambda=\sigma+i\omega$, where $\sigma\in\mathds{R}$ and $\omega\in\mathds{R}$, and computed the roots of the characteristic equations (D.6), (3.8) and (3.11) in the $(\sigma,\omega)$-plane using the Matlab subroutines roots and fsolve Mathworks .

We saw in Section 3.2 that the homeostasis steady states $\mathbf{X}^{*}_{2}$ and $\mathbf{Y}^{*}_{2}$ are locally asymptotically stable in all the versions of the Quartino model that we consider, and that the models are degenerate when $\beta=\gamma$, consequently here we will study the bifurcations that occur as parameters are varied from their homeostasis values with $\gamma>\beta$.

We begin by studying the Quartino model (2.2) starting from parameters used in Quartino2014 , so $n=4$, $a=k_{tr}$ and all the other parameters as in Table 1. We observed that changes to $\gamma$, the parameter relating the feedback of circulating G-CSF concentrations on the proliferating pool, and $a$, the transit rate between maturation compartments, can lead to a loss of stability giving rise to a periodic orbit via a Hopf bifurcation. Table 2 summarises the parameter pair values $(\gamma,a)$ necessary to induce such a loss in stability in $\mathbf{X}^{*}_{2}$, and the resulting period of the emerging periodic orbits.

If we let $\gamma$ be the bifurcation parameter and keep the remaining parameters at their homeostasis values (see Table 1), there is a Hopf bifurcation point at $\gamma=0.86766$. We verified numerically that $\mathbf{X}^{*}_{2}$ is locally asymptotically stable for $\beta<\gamma<0.86766$ and unstable if $\gamma>0.86766$. Of particular interest, as reflected in the bolded row of Table 2, we found a periodic orbit characteristic of cyclical neutropenia Dale2015 . Using the relation $a=n/$MMT with $n=4$ gives the value MMT$=123.33\,\text{hours}$, which is close to the mean maturation time of 133 hours for a patient under chemotherapy treatment reported by Quartino Quartino2014 .

In Figure 5 (left) we show the Hopf bifurcation curve for $\mathbf{X}^{*}_{2}$ on parameter space $(a,\gamma)$ for the Quartino model (2.2). The steady state is stable in the region below the Hopf curve and unstable otherwise. In Figure 5 (right) we see that increasing $\gamma$ for small values of the transit rate parameter $n$ lead to solutions with long period.

In the same vein, we also computed bifurcation points for the equilibrium $\mathbf{X}^{*}_{2}$ of distributed DDE model (2.31) using the characteristic equation (3.11) with $n=4$. As expected, given that the this model is simply a time-rescaling of the Quartino model (2.2) we obtain the same bifurcation points shown in Table 2.

In Figure 6 (left) we show the Hopf bifurcation curve for $\mathbf{Y}^{*}_{2}$ on parameter space $(n,\gamma)$ for the distributed DDE model (2.31). The steady state is stable in the region below the Hopf curve and unstable otherwise. Increasing $n$ lead to $\gamma$ and the period of the Hopf bifurcation converge to the bifurcation point of the discrete DDE model (2.33) shown in the bolded row of Table 4.

Table 4 reports bifurcation points computed for the steady state $\mathbf{Y}^{*}_{2}$ of the discrete DDE model (2.33) using the characteristic equation (3.8). Comparing the second rows of Tables 2 and 4, we note that the region of stability for the discrete DDE model $\beta<\gamma<0.69754$ is smaller than that of the distributed DDE with $n=4$, and of the equivalent ODE Quartino model $\beta<\gamma<0.86766$. Furthermore, we verified that in the limit $n\to\infty$ with $a=n/\tau$ or $a=(n+1)/\tau$ and holding $\tau$ fixed, the characteristic roots of the distributed DDE model (3.11) converge to the roots of the discrete DDE model (3.8) since (3.11) approaches to (3.8) when $n\to\infty$.

In Figure 7 (left) we show the Hopf bifurcation curve for $\mathbf{Y}^{*}_{2}$ on parameter space $(\tau,\gamma)$ for the discrete DDE model (2.33). The steady state is stable in the region below the Hopf curve and unstable otherwise. In Figure 7 (right) we see that there is an approximate linear relation between the period of the limit cycles and the mean value of the distributed delay $\tau$ for along all values of $\gamma$.

4.2 Bifurcations in the QSP model of granulopoiesis

We also studied whether parameter changes can lead the steady state $\mathbf{X}^{*}$ of system (2.34) to lose stability. We observed that changes in parameters related to proliferation and maturation lead to a loss of stability via a Hopf bifurcation, as reflected in Table 5. Inversely, we further verified the stability of the steady state when the half-maximal neutrophil proliferation constant satisfies $b_{N_{\!P}}\in[b_{N_{\!P}}\times{10}^{-3},b_{N_{\!P}}\times{10}^{3}]$, the rate of maturing neutrophil death satisfies $\gamma_{N_{\!M}}\in[\gamma_{N_{\!M}}\times{10}^{-2},\gamma_{N_{\!M}}\times{10}^{2}]$, and the neutrophil apoptosis rate in the bone marrow reservoir satisfies $\gamma_{N_{\!R}}\in[\gamma_{N_{\!R}}\times{10}^{-2},\gamma_{N_{\!R}}\times{10}^{2}]$.

We further investigated whether varying pairs of parameters in tandem could lead the steady state to lose stability via a Hopf bifurcation, as reflected in Table 6. An additional Hopf bifurcation point leading to an orbit of period of $20.94$ days was observed by changing four parameters simultaneously: $\eta_{N_{\!P}}^{h}=1.7\,\text{days}^{-1}$; $b_{N_{\!P}}=2.0$ ng/mL; $\eta_{N_{\!P}}^{\textit{min}}=1.3\,\text{days}^{-1}$; and $\tau_{N_{\!P}}=6.1\,\text{days}$.

5 The impact of stability and bifurcations on PK/PD considerations

In the PK/PD context, sensitivity analysis is frequently applied to investigate the impact of parameters variability on the system’s output. There the goal is to understand how predictions (outputs) change given changes to initial values (inputs). Still, when evaluating treatments, one may wonder how small changes to parameters affect the qualitative (e.g. existence and stability of equilibria, etc.) behaviour of the model. For example, if we have a priori information about a PK parameter’s variability and this parameter helps determine the model’s stability, bifurcation analysis can help to assess whether small changes within the range of the measured variability can bring about serious unintended shifts in the physiological system. Put another way, how are parameters changing when the system shifts stability or becomes unstable? Bifurcation analysis is rarely used in conventional PK/PD analyses, however the study of qualitative model behaviour is becoming increasingly recognised as an important tool for drug development Bakshi2016 ; Ghosh2013 . The potential impact of variability in PK parameters on the dynamics of the governing equations, which correspond to the PD aspects, is multifaceted. Within pharmacometrics, various situations must therefore be considered when assessing which (and how) parameters are susceptible to generating bifurcations when their values change.

In the simplest case, for physiological or drug parameters not influenced by drug concentration, no bifurcation can be generated through any PK variability. Examples of such parameters could include the maximal achievable response in an Emax model or a zero-order endogenous production rate. In that vein, we observe that in the Quartino model, the bifurcation point for $\gamma$ is not likely to be reached by realistic variations in G-CSF concentrations.

More familiar to the pharmacometrician is the case where different parameters values associated with specific cohorts and/or patient subpopulations correspond to individual states in the dynamical system. Here bifurcation analysis is analogous to Population-PK (Pop-PK) covariate analysis techniques used to separate and determine Pop-PK models for each subgroup. From this covariate analysis, one infers that for each state, there is a particular PK/PD model for which the (between subject) variability has been explained. Accordingly, Pop-PK covariate studies examine the impact of between subject variability is readily accounted for during the process of building a Pop-PK model. In the same way, bifurcation analysis of mathematical models helps to ascertain how changes to (variability in) model parameters affect the stable state. Of note, in the most complex case, where the parameters of a dynamical system are affected by PKs (which can also vary during therapy), within subject variability in PK is an important factor in determining model stability and a case-by-case model analysis must be carried out.

The impact of interindividual variability (IIV) is also an open question when considering the time-rescaled Quartino (2.30) and discrete delay (2.33) models. Since the discrete delay model is the equivalent to taking $n$ to infinity in (2.30), we investigated whether there would be altered behaviour due to IIV with increases in $n$. We began by generating 30 virtual patients using the docetaxel Pop-PK model of Bruno Bruno1996 and set Edrug=Slope$C_{\text{doc}}$, as in Quartino2014 , where $C_{\text{doc}}$ is the concentration of docetaxel in the central compartment. Using the same individual patient values for the generalised and time-rescaled Quartino models, we verified that solutions to (2.2) were identical as the solutions to (2.30) when the latter were rescaled according to (2.29) (not shown). We proceeded to compare the full Pop-PK predictions of (2.2) for these 30 virtual patients when $n=4$ and $n=10$. In Figure 8, increases to $n$ are not significantly affected by the inclusion of IIV, though increases to $n$ decidedly impact on the distribution of solutions. Given these results and to better visualise the impact of increasing $n$, we then compared the predictions of (2.2) to (2.33) (the case of infinite $n$) using only the typical parameter estimates. As seen in Figure 9, in the limit $n\to\infty$, solutions to the generalised Quartino model (and equivalently, the time-rescaled Quartino and distributed delay models) converge to that of the discrete delay. Further, as $n$ increases, so too does the dimension of the resulting ODE system for (2.2), which is not the case for the 3 equations of (2.33), thereby encouraging the use of a delay model to speed up simulation time. It should be noted that in all simulations described in this section, due to the misspecification of MMT as MMT=$(n+1)/k_{tr}$, to compare directly with the results of Quartino2014 $k_{tr}=k_{\textit{\tiny P}}$ were fixed to the value in Table 1, MMT was taken to be 133 hours, and $\tau$ was recalculate as $\tau=(4/5)\text{MMT}$, with $a=n/\tau$.

The analysis and results of Section 3.3 indicate the interest of performing stability and bifurcation studies in the pharmaceutical sciences setting. Through sensitivity analysis, we have previously concluded that models constructed from first-principles, such as (2.34), are robust to PKs by “sufficiently” accounting for the system’s physiological mechanisms Craig2016 . Thus average PK parameters without reference to the full Pop-PK model are reliably predictive of the behaviour of such physiological models. Here we extend our previous work to study how PK variability (variations in G-CSF concentrations) affect the stable states of the PD (physiological) system.

The granulopoietic system is indeed very robust around homeostatic parameter values. An increase in G-CSF concentrations corresponds to a decrease in the value of $\gamma_{Q}$, and, from Table 5, destabilisation of this equilibrium can only occur when $\gamma_{q}>\gamma^{h}_{Q}$, and thus such changes preserve homeostasis. With regard to changes in the proliferative progenitor compartment, either by variations in the parameters $\eta^{h}_{N_{P}}$ or $b_{v}$, physiologically realistic scenarios preclude reaching the bifurcation values: for example, in the case of $\eta^{h}_{N_{P}}$, artificial removal of G-CSF from the body would be required.

6 Discussion

Mathematical pharmacology, defined as the study of mathematical approaches to pharmacological processes, is increasingly recognised as a quantitative methodology critical to understanding pharmaceutical treatments and their efficacy while simultaneously raising compelling mathematical problems vanderGraaf2016 . Using granulopoiesis as a backdrop, in the present paper we have examined the connections between the familiar PK/PD model formalism originally proposed by Friberg Friberg2002 and adapted by Quartino Quartino2014 and a discrete delay model of neutrophil production, connected via a distributed delay model. Crucially, we have shown how the stability of each model can be studied straightforwardly via this latter distributed delay model, underlining the advantage of being able to transfer between these equivalent expressions and motivating the present analysis. We examined the impact of the inclusion of IIV on the solutions to each of the models, and determined that variations are driven through increases to $n$ rather than the presence of variability; as $n\to\infty$, solutions of the generalised (time-rescaled) Quartino transit compartment model converge to that of the time-rescaled discrete DDE model. Last, using our previously published QSP model of the negative feedback relationship between granulopoiesis and G-CSF, we have identified several Hopf bifurcations through bifurcation analysis, a technique not commonly applied in the classical PK/PD analyses, and reviewed the impact the interpretation of such bifurcations can have on our understanding of pharmacological systems when used in concert with more common sensitivity and variability analyses. We would like to highlight two results in particular. First, the distributed delay model (2.31) exhibits wider regions of stability around the steady state $\mathbf{Y}_{2}^{*}$ as compared to the discrete DDE model (2.33), consistent with the result that “distributed delays are inherently more stable than the same system with discrete delays” Campbell2009 . Second, we identified Hopf bifurcations in the distributed and discrete delay forms of the Quartino model with periods corresponding to those in cyclic neutropenic patients by varying the feedback parameter $\gamma$ and the delay $\tau$, demonstrating how bifurcation analyses can be applied in mathematical pharmacology to understand the pathogenesis towards diseases.

Perhaps the most immediately consequential conclusion drawn here is the incorrect definition of the mean transit/maturation time in the original and subsequent applications of the Friberg model. As previously mentioned, by setting $k_{\textit{\tiny P}}=a=k_{tr}$ with $G(t)=G_{0}$, the MTT(MMT) was originally expressed as $(n+1)/k_{tr}$. However, we have shown that the mean delay of the distributed delay model is constant and instead given by $\tau=n/a$ (which, when $a=k_{tr}$ is then clearly given by $n/k_{tr}$). Thus it is mathematically incorrect to set MMT$=(n+1)/a$ as it treats the proliferative pool as an additional transit compartment and this formulation cannot be recovered via the linear chain technique. The generalised Quartino model (2.2) explicitly decouples the maturation time and the production rate of cells to eliminate this problem. Additionally, we highlighted the mathematical issue presented when $a$ is non-constant, as in Quartino2014 , to the derivation (2.14), which is essential to the linear chain technique to recover the correct ODE formulation from the distributed delay model (2.26).

Further, since, in Friberg2002 and its various extensions and applications, the parameter $k_{tr}$ is determined via the MMT, the mean maturation time is fit and then the rate of transit through each compartment is determined via the equation MMT$=(n+1)/k_{tr}$. This leads to disparate estimates for the maturation process, ranging, for example, from 102 hours ($n=6$) in Quartino2012 to 210 hours ($n=4$) in Quartino2014 . Physiological labelling studies report a much narrower range of maturation times (6.4 days in Price and 6.9 days in Dancey1976 , for example). Thus, allowing the MMT to vary widely is not physiologically consistent and further introduces additional mathematical difficulties since, in general, $\tau=n/a$, where $a$ is not necessarily equal to $k_{tr}$ nor $k_{\textit{\tiny P}}$.

We therefore emphasise that this work provides further motivation to systematically incorporate, from first principles, the physiological architecture yielding the proper mathematical formulation of pharmacological models.

Acknowledgements

DCS was supported by National Council for Scientific and Technological Development of Brazil (CNPq) postdoctoral fellowship 201105/2014-4. MC was supported by an Natural Sciences and Engineering Research Council of Canada (NSERC) postdoctoral fellowship and grant DP5OD019851 from the Office of the Director at the National Institutes of Health to her PI. TC was supported by the Alberta government via the Sir James Lougheed Award of Distinction as well as the Centre de Recherche Mathématiques, Montréal. FN and JL are funded by FN’s NSERC Industrial Chair in Pharmacometrics, supported by Novartis, Pfizer, and inVentiv Health Clinics, and an FQRNT projet d’équipe. JB and ARH are grateful to NSERC for funding through the Discovery Grant program. We are appreciative for our many very useful discussions with Michael C. Mackey.

Appendix A Time Rescaling of Quartino Model

Here we show how the time rescaling (2.29) we applied to the generalised Quartino model (2.2) relates to the state-dependent delays that are used in the QSP granulopoiesis model (2.34).

For the Friberg model (2.1) we have average maturation delay $\tau$ given by $\tau=n/a$. Hence

[TABLE]

For the Quartino model (2.2) the maturation rate $a$ is replaced by $a(G/G_{0})^{\beta}$ and hence the time-dependent maturation delay $\alpha(t)$ for this model is given by

[TABLE]

Thus the mean maturation time $\tau$ for the time-rescaled Quartino model (2.30) is related to $\alpha(t)$ by

[TABLE]

Equation (A.1) defines $\alpha(t)$ by a threshold condition. This is completely analogous to the threshold condition (2.38) used to define the state-dependent delay $\tau_{N_{\!M\!}}(t)$ in the QSP model (2.34).

Differentiating (A.1) using Leibniz rule we obtain an expression for the evolution of $\alpha(t)$ as

[TABLE]

which can be rewritten as

[TABLE]

and determines the evolution of $\alpha(t)$. An analogous expression was derived in Craig Craig2016c for the evolution of $\tau_{N_{\!M\!}}(t)$ in the QSP model (2.34).

Appendix B Positivity of Solutions

We show positivity of the solutions to the models considered in the paper.

B.1 Positivity of solutions to the Quartino endogenous G-CSF model

Consider the generalised Quartino model (2.2).

Lemma 1

Assume that the parameters in (2.2) are strictly positive and that $G(0)\geqslant 0$.

1. i)

Then $G(t)>0$ for all time $t>0$. 2. ii)

If $N(t)\geqslant 0$ for $t\in[0,t_{\textit{\tiny N}}]$ then $G(t)\leqslant\max\{k_{in}/k_{e},G(0)\}$ for all $t\in[0,t_{\textit{\tiny N}}]$.

Proof

i) Either $G(0)>0$ or $G(0)=0$ and $\frac{\textrm{d}G}{\textrm{d}t}(0)=k_{in}>0$. In both cases $G(t)>0$ for $t\in(0,\varepsilon)$ for some $\varepsilon>0$. Assume, for contradiction, that there exists $t_{\textit{\tiny G}}>0$ such that $G(t_{\textit{\tiny G}})=0$ but $G(t)>0$ for $t\in(0,t_{\textit{\tiny G}})$. Since $G$ is decreasing at $t=t_{\textit{\tiny G}}$ this implies that $\frac{\textrm{d}G}{\textrm{d}t}(t_{\textit{\tiny G}})\leqslant 0$. But this is contradicted by (2.2) which implies that $\frac{\textrm{d}G}{\textrm{d}t}(t_{\textit{\tiny G}})=k_{in}>0$ if $G(t_{\textit{\tiny G}})=0$. Thus there exists no such that time $t_{\textit{\tiny G}}$, and hence $G(t)>0$ for all $t>0$. ∎

ii) For $t\in[0,t_{\textit{\tiny N}}]$ we have $\frac{\textrm{d}G}{\textrm{d}t}\leqslant k_{in}-k_{e}G$, and hence $\frac{\textrm{d}G}{\textrm{d}t}<0$ if $G(t)>k_{in}/k_{e}$. The result follows. ∎

Lemma 2

Assume that the parameters in (2.2) are strictly positive, and that the initial conditions satisfy $P(0)>0$ and $G(0)\geqslant 0$. Furthermore, assume that there exists $t_{\textit{\tiny N}}>0$ such that $N(t)\geqslant 0$ for $t\in[0,t_{\textit{\tiny N}}]$. Then $P(t)>0$ for all time $t\in[0,t_{\textit{\tiny N}}]$.

Proof

By Lemma 1 for $t\in[0,t_{\textit{\tiny N}}]$ we have $0\leqslant G(t)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$. Thus from (2.2) we have

[TABLE]

which (using the continuous Gronwall lemma) implies that

[TABLE]

Lemma 3

Assume that the parameters in (2.2) are strictly positive, and that the initial conditions satisfy $P(0)>0$, $G(0)\geqslant 0$ and $T_{j}(0)\geqslant 0$ for $j=1,2,\ldots,n$. Then there exists $\varepsilon>0$ such that $T_{j}(t)>0$ for all time $t\in(0,\varepsilon)$ for all $j=1,2,\ldots,n$. Furthermore, if there exists $t_{\textit{\tiny N}}>0$ such that $N(t)\geqslant 0$ for $t\in[0,t_{\textit{\tiny N}}]$ then $T_{j}(t)>0$ for all time $t\in(0,t_{\textit{\tiny N}}]$ for all $j=1,2,\ldots,n$.

Proof

We proceed by induction on $j$. For $j=1$, either (i) $T_{1}(0)>0$ and $G(0)\geqslant 0$, or (ii) $T_{1}(0)=0$ and $G(0)>0$, or (iii) $T_{1}(0)=G(0)=0$. In case (ii) $\frac{\textrm{d}T_{1}}{\textrm{d}t}(0)=k_{tr}\left(\frac{G(0)}{G_{0}}\right)^{\beta}P(0)>0$, while in case (iii) $T_{1}(0)=\frac{\textrm{d}T_{1}}{\textrm{d}t}(0)=0$.

For all three cases, $P(0)>0$ implies that that there exists $\varepsilon>0$ such that $P(t)>0$ for $t\in(0,\varepsilon)$. Lemma 1 ensures also that $G(t)>0$ for $t\in(0,\varepsilon)$. But now, for $t\in(0,\varepsilon)$ if $T_{1}(t)\leqslant 0$ we have the strict inequality

[TABLE]

Thus if there exists $t_{1}\in(0,\varepsilon)$ such that $T_{1}(t_{1})\leqslant 0$ then $\frac{\textrm{d}}{\textrm{d}t}T_{1}(t)>0$ for all $t<t_{1}$ and therefore $T_{1}(0)<0$, which contradicts the initial condition. Hence there exists no such time $t_{1}$, and so $T_{1}(t)>0$ for $t\in(0,\varepsilon)$.

For general $j$, it is shown that $T_{j}(t)>0$ for $t\in(0,\varepsilon)$ similarly. The positivity of $G(t)>0$ and $T_{j-1}(t)>0$ for $t\in(0,\varepsilon)$ along with the strict inequality

[TABLE]

for all $t\in(0,\varepsilon)$ if $T_{j}(t)\leqslant 0$ similarly ensures that actually $T_{j}(t)>0$ for $t\in(0,\varepsilon)$, which establishes the result for any finite $n$.

Now, from Lemmas 1 and 2, if there exists $t_{\textit{\tiny N}}>0$ such that $N(t)\geqslant 0$ for $t\in[0,t_{\textit{\tiny N}}]$ then $P(t)>0$ for all time $t\in[0,t_{\textit{\tiny N}}]$ and $0\leqslant G(t)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$. Now we have $T_{1}(t)>0$ for $t\in(0,\varepsilon)$ and for $j=1$

[TABLE]

which, similar to the proof of Lemma 2, implies $T_{1}>0$ for all $t\in(0,t_{\textit{\tiny N}}]$. But if $T_{j-1}(t)>0$ for all $t\in(0,t_{\textit{\tiny N}}]$ then $T_{j}(t)>0$ satisfies (B.1) and it follows by induction that $T_{j}(t)>0$ for all $t\in(0,t_{\textit{\tiny N}}]$ for $j=1,\ldots,n$ if $N(t)\geqslant 0$ for $t\in[0,t_{\textit{\tiny N}}]$. ∎

Theorem B.1

Assume that the parameters in (2.2) are strictly positive, and that the initial conditions satisfy $P(0)>0$, $G(0)\geqslant 0$, $N(0)\geqslant 0$ and $T_{j}(0)\geqslant 0$ for $j=1,2,\ldots,n$. Then $P(t)>0$, $0<G(t)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$, $N(t)>0$ and $T_{j}(t)>0$ for $j=1,2,\ldots,n$ for all $t>0$.

Proof

Lemma 1 implies that $G(t)>0$ for all $t>0$. While $P(0)>0$ and Lemma 3 imply that $P(t)>0$ and $T_{j}(t)>0$ for $j=1,2,\ldots,n$ for $t\in(0,\varepsilon)$.

Now similar to the proof of Lemma 3,

[TABLE]

implies that $N(t)>0$ for $t\in(0,\varepsilon)$. Then while $N(t)>0$, we have $G(t)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$ and

[TABLE]

But this last inequality implies that $N(t)>0$ for all $t>0$. It follows from Lemmas 1, 2 and 3 that $G(t)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$, $P(t)>0$, and $T_{j}(t)>0$ for $j=1,\ldots,n$ for all $t>0$. ∎

Although an essential part of the solution positivity proofs in this section was to show that solutions decay with at most a bounded linear rate for general parameters, these are not the dynamics that we actually expect to observe. For normal individuals/subjects we should have $\gamma>\beta$ (as is the case for the parameters in Table 1), in which case when $G(t)$ is sufficiently large $\frac{\textrm{d}P}{\textrm{d}t}$ is positive, as is required for the feedback loops to function effectively.

Since the Quartino model (2.2) is equivalent to the time-rescaled Quartino model (2.30) and the distributed DDE (2.31), positivity of solutions to those models follows directly from Theorem B.1 (it is important to note that $\hat{t}(t)>0$ for all $t>0$ follows from (2.29) and the positivity of $G(t)$). However, this only establishes positivity of solutions for the distributed DDE (2.31) when $n$ is an integer. When can show positivity directly for both the distributed DDE (2.31) for general real $n$ and also for the discrete DDE (2.33).

Theorem B.2

Assume that the parameters in distributed delay DDE (2.31) or the discrete delay DDE (2.33) are strictly positive, and that the initial conditions satisfy $G(0)>0$, $N(0)\geqslant 0$, and $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)=\varphi(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)$ for $\hat{t}\leqslant 0$ where $\varphi$ is continuous and $P(0)=\varphi(0)>0$. Then Then $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)>0$ and $0<G(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)\leqslant M=\max\{k_{in}/k_{e},G(0)\}$ for all $\hat{t}>0$. Finally $N(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)>0$ for all $\hat{t}>0$ for the system (2.31), and for all $\hat{t}\>\tau$ for the system (2.33).

Proof

The proof uses similar ideas to the proof of Theorem B.1, so we just outline the details here. We have $P(0)>0$, and similar to the proof of Lemma 2 the rate of decrease of $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)$ is bounded so $P(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)>0$ for all $\hat{t}>0$. But now for (2.31), the positivity of $P$ ensures that $\frac{\textrm{d}}{\textrm{d}t}N>0$ if $N=0$, which ensures that $N(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)>0$ for all $\hat{t}>0$. For the model (2.33), $P(\hat{t}-\tau)>0$ for $\hat{t}\geqslant\tau$ leads to the positivity of $N(\hskip 0.50003pt\hat{t}\hskip 0.50003pt)>0$ for all $\hat{t}>\tau$. Finally the bounds on $G$ are derived similarly to Lemma 1. ∎

B.2 Positivity of solutions to the QSP granulopoiesis model

Consider the QSP granulopoiesis model (2.34). Using the constraints listed in (C.1) and setting

[TABLE]

yields

[TABLE]

Lemma 4

Consider the initial value problem $\mathcal{P}$, given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35), with $G_{1,0}+G_{2,0}>0$. Assume that $G_{\!\textit{prod}},k_{ren}$ and $k_{int}$ are positive constants and set $\beta=\max\{k_{ren},k_{int}\}$. Then

[TABLE]

Proof

A simple calculation shows

[TABLE]

Multiplying by $e^{\beta t}$ and rearranging gives

[TABLE]

Integrating the inequality from [math] to $t$ yields

[TABLE]

Taking $K=[G_{1,0}+G_{2,0}]$ yields the claim. ∎

Rearranging the bound of Lemma 4 gives:

[TABLE]

Lemma 5

Consider the initial value problem, $\mathcal{P}$, given equations (2.34a)– (2.34f) and the ICs and histories given in (2.35), with $G_{1,0}+G_{2,0}>0$ and $G_{1,0}\geqslant 0$. Assume that $G_{\!\textit{prod}},k_{ren},k_{21}$ and $k_{int}$ are strictly positive constants. Then, $G_{1}(t)>0$ for all $t>0$.

Proof

Assume, for contradiction, that there exists a time $s>0$ such that $G_{1}(s)<0$. As $G_{1}$ is a solution of the differential equation, it is continuously differentiable. By the IVT, there must exist a time $T^{*}$ such that $G_{1}(T^{*})=0$ with $s>T^{*}$. Using (B.3), at $t=T^{*}$:

[TABLE]

The Mean Value Theorem (MVT) yields a contradiction to $G_{1}(s)<0$. ∎

Lemma 6

Consider the initial value problem $\mathcal{P}$, equations (2.34a)– (2.34f) and the ICs and histories given in (2.35), with $G_{1,0}+G_{2,0}>0$, $G_{1,0}\geqslant 0$ and $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau_{Q},0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$. Assume that $G_{\!\textit{prod}},k_{ren},k_{21}$ and $k_{int}$ are strictly positive constants. Then, $Q(t)>0$ for all time $t>0$.

Proof

Assume, for contradiction, that there exists a $s^{*}$ such that $Q(s^{*})=0$. The assumption on $\varphi_{1}$ and the continuity of $Q$ allows us to assume that $Q(s)>0$ for all $s<s^{*}$.

The constraints in (B.2) ensures that $\beta(Q(t-\tau_{Q}))>0$.

Then:

[TABLE]

Using the result of Lemma 5 and that $\kappa(G_{1}(t))$ is a monotonically increasing function for $G_{1}(t)>0$ gives

[TABLE]

Setting $\xi=\kappa^{max}+\kappa_{\delta}+f_{Q}$ gives:

[TABLE]

An argument similar to that in Lemma 4 ensures that $Q(t)>0$ for all $t>0$. This contradicts the existence of $s^{*}$ and establishes the positivity of $Q(t)$ for all $t>0$. ∎

Lemma 7

Consider the initial value problem $\mathcal{P}$, given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35), with $G_{1,0}+G_{2,0}>0$, $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and $N_{0}+N_{R,0}>0$. Assume that $G_{\!\textit{prod}},k_{ren},k_{21}$ and $k_{int}$ are strictly positive constants. Then $N_{\!R}(t)+N(t)>0$ for all $t>0$.

Proof

A simple calculation gives:

[TABLE]

Applying Lemmas 5, 6 and (B.2), the following holds

[TABLE]

Define $\alpha=\textrm{max}[\gamma_{N_{\!R}},\gamma_{N}]$ and calculate:

[TABLE]

An argument similar to that of Lemma 4 gives the positivity of the sum. ∎

Lemma 8

Consider the initial value problem $\mathcal{P}$ given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35) such that $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and

[TABLE]

Assume that $k_{21}$ and $k_{int}$ are strictly positive constants. Then, $G_{2}(t)\leqslant V[N_{\!R}(t)+N(t)]$ for all $t>0$.

Proof

Assume, for contradiction, that there is a time $s>0$ with $G_{2}(s)>V[N_{\!R}(t)+N(t)]$. Then, there must exist a time $T^{*}$ such that $G_{2}(T^{*})=V[N_{\!R}(T^{*})+N(T^{*})]$. At $t=T^{*}$:

[TABLE]

Lemma 7 ensures that $G_{2}(T^{*})=V[N_{\!R}(T^{*})+N(T^{*})]>0$. Therefore $\frac{\textrm{d}}{\textrm{d}t}G_{2}(t)|_{t=T^{*}}<0$, which contradicts the MVT and there can be no $s$. ∎

Lemma 9

Consider the initial value problem $\mathcal{P}$ given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35) such that $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and

[TABLE]

Assume that $G_{\!\textit{prod}},k_{ren},k_{21}$ and $k_{int}$ are strictly positive constants. Then, $G_{2}(t)>0$ for all $t>0$.

Proof

Using Lemma 5 and the bound from Lemma 8 gives

[TABLE]

and

[TABLE]

An argument similar to Lemma 4 gives the result. ∎

Lemma 10

Consider the initial value problem $\mathcal{P}$ given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35) such that $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and

[TABLE]

Assume that $k_{21}$ and $k_{int}$ are positive constants. Then, $N_{\!R}(t)>0$ for all $t>0$.

Proof

The function $A_{N}(t)$ is positive for all time, the constraints in (B.2), the assumptions on $\varphi_{1,2}$ and Lemmas 6 and 5 ensure that $A_{N}(t)\kappa(G_{1}(t-\tau_{N}(t)))Q(t-\tau_{N}(t))\frac{V_{\!N_{\!M\!}}(G_{1}(t))}{V_{\!N_{\!M\!}}(G_{1}(t-\tau_{N_{\!M\!}}(t)))}\geqslant 0$. Then

[TABLE]

Lemma 8 and the positivity of $N(t)+N_{\!R}(t)$ and $G_{2}(t)$ bounds $G_{BF}\in[0,1]$. Therefore, $\varphi_{N_{\!R}}(G_{BF})$ is a continuous function on a compact domain and is therefore bounded below by 0 and above by $\varphi_{N_{\!R}}^{max}$. Then

[TABLE]

and an argument similar to Lemma 4 yields the result. ∎

Lemma 11

Consider the initial value problem $\mathcal{P}$ given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35) such that $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and

[TABLE]

Assume that $k_{21}$ and $k_{int}$ are strictly positive constants. Then, $N(t)>0$ for all $t<0$.

Proof

Using (B.2) and the Lemmas 5 and 10 to ensure the positivity of $N_{\!R}(t)$ and $\varphi_{N_{\!R}}$ yields:

[TABLE]

A similar argument to that used in Lemma 4 yields the result. ∎

Together, these results lead to the following theorem.

Theorem B.3

Consider the initial value problem $\mathcal{P}$ given by equations (2.34a)– (2.34f) and the ICs and histories given in (2.35) such that $\varphi_{1}(s)\geqslant 0$ for $s\in[-\tau,0]$ with $\varphi_{1}(s^{*})>0$ for at least one $s^{*}$, $\varphi_{2}\geqslant 0$ for $s\in[-\tau,0]$ and

[TABLE]

with $I_{G}(t)=0$. Moreover, assume that (C.1) is satisfied along with strictly positive model parameters. Finally, assume that the history functions $\varphi_{1,2}$ are positive at least once and are non-negative in their domain. Then

the solutions of $\mathcal{P}$, $\mathbf{x}(t)=\left(Q(t),N_{\!R}(t),N(t),G_{1}(t),G_{2}(t)\right)$ remain component wise positive for all time. 2. 2.

the solutions of $\mathcal{P}$, $\mathbf{x}(t)=\left(Q(t),N_{\!R}(t),N(t),G_{1}(t),G_{2}(t)\right)$ remain component wise bounded for all time.

Proof

1. The results of Lemmas 5; 6; 9; 10; and 11 give the result.

2. In Mackey1994 , the authors prove that solutions of the equation:

[TABLE]

are bounded above by a finite $x_{1}$. Setting $\tilde{\delta}=\kappa^{min}+\kappa_{\delta}$ gives:

[TABLE]

which implies that $Q(t)$ is bounded by some $Q_{1}$.

$V_{N_{M}}$ is an increasing function bounded above by $V_{max}$ and bounded below by $V(0)>0$. Finally,

[TABLE]

This gives

[TABLE]

Then, setting $K=\exp[\eta_{N_{\!P}}^{max}\tau_{N_{\!P}}]f_{0}Q_{1}\frac{V_{max}}{V_{N_{M}}(0)}$ and $\beta=\textrm{min}[\gamma_{N},\gamma_{N_{\!R}}]$ gives

[TABLE]

Equation (B.5) is equivalent to

[TABLE]

Integrating this inequality gives

[TABLE]

As $N_{R},N$ are both positive, they are individually bounded by this constant.

Now,

[TABLE]

Setting $\alpha=\textrm{min}[k_{int},k_{ren}]$ gives:

[TABLE]

A similar argument to above gives

[TABLE]

Again, as both $G_{1}(t)$ and $G_{2}(t)$ are positive, they are individually bounded by this constant.

Then, solutions of the mathematical model $\mathcal{P}$ with positive initial conditions remain positive and bounded. ∎

Appendix C Existence and uniqueness of the Homeostatic Steady State of the QSP granulopoiesis model

Consider the QSP granulopoiesis model defined by the system of DDEs (2.34), and associated initial conditions and histories, with strictly positive parameters satisfying the constraints

[TABLE]

The following results will demonstrate the existence and uniqueness of a homeostasic steady state of this model.

Proposition 1

Assume that

[TABLE]

and define

[TABLE]

If $\overline{G}_{1}>G_{\!\textit{prod}}$, then the system of DDEs given by (2.35) has a unique positive homeostatic steady state.

Proof

Consider the differential equation given in (2.35). Define

[TABLE]

The steady states are given implicitly by the solutions of

[TABLE]

and are denoted by $Q^{h},N_{\!R}^{h},N^{h},G_{1}^{h},G_{2}^{h}$ respectively. Moreover, denote the homeostatic values of each function with the superscript $h$.

A simple calculation to find the non-zero solution of

[TABLE]

gives the expression for the homeostatic steady state of $Q(t)$

[TABLE]

The conditions in the statement of Proposition 1 ensure that the steady state is positive. Using $Q^{h}$, we calculate the solution $N_{\!R}^{h}$ of

[TABLE]

The homeostatic value is

[TABLE]

The final homeostasis value that can be easily expressed is the solution of

[TABLE]

Using the homeostatic value of $N_{\!R}^{h}$ and solving the steady state equation gives

[TABLE]

Now, consider the coupled G-CSF kinetics given by equations (2.34e) and (2.34f) whose homeostatic sum is expressed as

[TABLE]

Isolating $G_{2}^{h}$ as a function of $G_{1}^{h}$ and substituting into (2.34f)

[TABLE]

gives the following exponential polynomial in $G_{1}^{h}$.

[TABLE]

where

[TABLE]

Then, any root of $F$ and the corresponding $G_{2}$ value from (C.2) will be a solution of (C.3). The condition

[TABLE]

guarantees that the coefficients are positive.

A simple calculation shows that

[TABLE]

By the IVT, there must exist a solution $G_{1}^{h}\in(0,\zeta/\delta)$ such that

[TABLE]

Since for all $x\in\mathds{R}^{+}$

[TABLE]

thus $F$ is a strictly increasing function and the zero found in (C.5) is unique.

Therefore, the corresponding $G_{2}^{h}$ given by (C.2) is unique. The condition $0<G_{1}^{h}<\zeta/\delta$ together with (C.2) and $\zeta/\delta=G_{\!\textit{prod}}/k_{ren}$ ensures the positivity of $G_{2}^{h}$. Thus the homeostatic steady state $(Q^{h},N_{R}^{h},N^{h},G_{1}^{h},G_{2}^{h})$ exists, is positive and unique. ∎

Appendix D Stability of the equilibria

D.1 Characteristic equations of the Quartino model

The characteristic equations for the time-rescaled Quartino model (2.30), for the equivalent distributed delay DDE (2.31), and for the discrete delay DDE model (2.33) have relatively simple form and are stated explicitly in the main text as equations (3.10) and (3.8). In contrast the Quartino model (2.2) has a much more complicated form both for general $n$ and when $n=4$ as in Quartino Quartino2014 , which we derive here.

The Quartino model (2.2) has characteristic equation given by (3.11) where $\mathds{J}$ is the Jacobian of $\mathbf{F}$ in (3.2). Differentiating the terms in (2.2) for general $\mathbf{X}$ we obtain the Jacobian

[TABLE]

where

[TABLE]

Then the characteristic function for the Quartino model (2.2) is the $(n+3)$-degree polynomial in $\lambda$:

[TABLE]

This is considerably more complicated than the characteristic function of the time-rescaled Quartino model which we also stated for general $\mathbf{X}$ in equation (3.11). Fortunately (D.2) does simplify somewhat at the steady states. At both $\mathbf{X}_{1}^{*}$ and $\mathbf{X}_{2}^{*}$ we have

[TABLE]

and hence

[TABLE]

Additionally at $\mathbf{X}_{2}^{*}$ only, $G=G_{0}$ implies $J_{11}=1$, while at $\mathbf{X}_{1}^{*}$ only, $P=T_{n}=0$ implies that $T_{1,n+3}=T_{n+2,n+3}=0$. So

[TABLE]

with

[TABLE]

and

[TABLE]

At $\mathbf{X}_{2}^{*}$ from (D.2) we obtain

[TABLE]

Equations (D.5) and (D.6) give the characteristic function as a polynomial of degree $n+3$ at the two steady states $\mathbf{X}_{1}^{*}$ and $\mathbf{X}_{2}^{*}$ of the generalised Quartino model (2.2) for general values of the parameters, including for the parameters used in Quartino Quartino2014 when $n=4$ and the polynomials have degree $7$.

D.2 Stability of the equilibria of the Quartino ODE model

The stability of the steady states of the Quartino model (2.2) was already considered in Theorem 3.1, which was proved using the alternate time-rescaled forms of the model. Here we present proofs of stability results directly in the original formulation, mainly to demonstrate how much more difficult they are.

Proposition 2

The equilibrium point $\mathbf{X}^{*}_{1}$ of the Quartino model (2.2) is locally asymptotically stable when $\gamma<\beta$ and unstable when $\gamma>\beta$.

Proof

For the case $\beta>\gamma$, we show asymptotic stability of $\mathbf{X}^{*}_{1}$ using a Lyapunov function Meiss2007 applied to the first $(n+2)$ coordinates of the system.

First consider the dynamics of $G(t)$. Let $\lambda\in(0,1)$ and suppose $N(t)\in[0,\lambda N_{0}]$ for $t\in[0,t_{\textit{\tiny N}}]$, and let

[TABLE]

Then from (2.2e) we have $\frac{\textrm{d}}{\textrm{d}t}G\geqslant 0$ when $G/G_{0}<G_{\textit{ratio}}$, and the positivity of $N(t)$ implies that $\frac{\textrm{d}}{\textrm{d}t}G\leqslant 0$ when $G>k_{in}/k_{e}$. Hence choosing $G(0)\in[G_{0}G_{\textit{ratio}},(1+\lambda)k_{in}/k_{e}]$, ensures that $G(t)\in[G_{0}G_{\textit{ratio}},(1+\lambda)k_{in}/k_{e}]$ for all $t\in[0,t_{\textit{\tiny N}}]$, under the assumption that $N(t)\in[0,\lambda N_{0}]$ for $t\in[0,t_{\textit{\tiny N}}]$. Now let

[TABLE]

where the parameters satisfy $\mu>0$ and $\alpha>0$, with values to be specified below. Note that $V(\mathbf{X}^{*}_{1})=0$, while $V(\mathbf{X})>0$ if any of $P$, $T_{j}$, $N$ is non-zero. Then differentiating and using (2.2) we obtain

[TABLE]

Notice that all the terms on the right-hand side are non-positive, except possibly for the first and last term. But with $k_{tr}=k_{\textit{\tiny P}}$, $\beta>\gamma$ and $G(t)/G_{0}\geqslant G_{\textit{ratio}}>1$ it follows that

[TABLE]

and hence for $\mu>0$ sufficiently large the coefficient of $P(t)^{2}$ is strictly negative. Similarly, since $G(t)$ is bounded above for $t\in[0,t_{\textit{\tiny N}}]$ for $\alpha>0$ sufficiently small the coefficient of $N(t)^{2}$ is also strictly negative. It follows that $\frac{\textrm{d}}{\textrm{d}t}V(\mathbf{X}(t))<0$ unless $P=T_{j}=N=0$.

If $V(\mathbf{X}(t))<\frac{\alpha\lambda^{2}N_{0}^{2}}{2}$, then (D.7) implies that $N(t)<\lambda N_{0}$. Therefore, we choose initial conditions such that $G(0)\in[G_{0}G_{\textit{ratio}},(1+\lambda)k_{in}/k_{e}]$ and

[TABLE]

Then since $V(\mathbf{X}(t))$ is nonincreasing, we have $N(t)<\lambda N_{0}$ for all $t>0$, then $G(t)\in[G_{0}G_{\textit{ratio}},(1+\lambda)k_{in}/k_{e}]$ for all $t>0$ and $V(\mathbf{X}(t))\to 0$ as $t\to\infty$ which implies that $[P(t),T_{1}(t),\ldots,T_{n}(t),N(t)]\to[0,0,\ldots,0,0]$ as $t\to\infty$.

Finally considering (2.2e), we have $\frac{\textrm{d}}{\textrm{d}t}G<0$ if $G>k_{in}/k_{e}$ and $\frac{\textrm{d}}{\textrm{d}t}G>0$ if $G<k_{in}/(k_{e}+k_{\textit{\tiny ANC}}N(t))$ with $N(t)\to 0$ as $t\to\infty$, hence $G(t)\to k_{in}/k_{e}$ as $t\to\infty$. This completes the proof of local asymptotic stability of $\mathbf{X}^{*}_{1}$ when $\gamma<\beta$.

For the case $\gamma>\beta$, we show that $\mathbf{X}^{*}_{1}$ is unstable, using linearization theory, by showing there is a real positive characteristic root in this case. The characteristic function evaluated at $\mathbf{X}^{*}_{1}$ is given by (D.5) and statisfies $\Delta_{1}(\lambda)\to+\infty$ as $\lambda\to+\infty$. Hence, by the IVT, to show that there exists a positive eigenvalue, $\lambda>0$ such that $\Delta_{1}(\lambda)=0$, it is sufficient to show that $\Delta_{1}(0)<0$. But from (D.5)

[TABLE]

From the constraint in equation (2.8), we have that $k_{in}=G_{0}(k_{e}+k_{\textit{\tiny ANC}}N_{0})>G_{0}k_{e}$, then $\gamma>\beta$ and the positivity of all the parameters imply that $\Delta_{1}(0)<0$ as required. Thus $\Delta_{1}(\mathbf{X}^{*}_{1})$ has a characteristic value $\lambda>0$ and $\mathbf{X}^{*}_{1}$ is unstable. ∎

Proposition 3

The equilibrium point $\mathbf{X}^{*}_{2}$ of the Quartino model (2.2) is unstable when $\gamma<\beta$.

Proof

Consider the characteristic function $\Delta_{2}(\lambda)$ evaluated at $\mathbf{X}_{2}^{*}$. From (D.6) we have $\Delta_{2}(\lambda)\to+\infty$ as $\lambda\to\infty$, while recalling that $J_{11}=0$ and using (2.7) and (2.8)

[TABLE]

From (2.8) we have $k_{e}<(k_{e}+k_{\textit{\tiny ANC}}N_{0})=k_{in}/G_{0}$, thus $\Delta_{2}(0)<0$. Thus the IVT guarantees that (D.6) has at least one real positive root and therefore $\mathbf{X}^{*}_{2}$ is unstable. ∎

D.3 Characteristic Equation of the discrete delay DDE Quartino model

The discrete delay DDE model (2.33), obtained as the $n\to\infty$ limit of the time-rescaled Quartino model (2.30) has characteristic equation (3.7). The $3\times 3$ linearisation matrices $\mathds{A}$ and $\mathds{B}$ from (3.6), are derived by differentiating the terms on the right-hand side of (2.33). We will denote their entries by $A_{ij}$ and $B_{ij}$ for $i,j=\{1,2,3\}$. Let us define the ratio between G-CSF concentrations $G_{r}\coloneqq G^{*}/G_{0}$ and let $\mathbf{Y}^{*}=(P^{*},N^{*},G^{*})$ be a generic steady state then

[TABLE]

where

[TABLE]

Using the matrices (D.8) we rearrange the determinant from equation (3.7) and obtain

[TABLE]

from where we get the characteristic equation (3.8) with the coefficients $a_{2}$, $a_{1}$, $a_{0}$ and $b$ given by

[TABLE]

Recalling the steady states $\mathbf{Y}^{*}_{1}$ and $\mathbf{Y}^{*}_{2}$ given, respectively, by (3.16) and (3.17), calculated with the parameter constraints (2.4) and (2.8), we can evaluate the coefficients in (D.9) at each steady state. For $\mathbf{Y}^{*}_{1}$ we have $P^{*}=N^{*}=0$ which implies that $A_{13}=A_{23}=0$ and we have

[TABLE]

For $\gamma>\beta$ we have $a_{0}<0$ in (D.10), since $k_{in}>k_{e}G_{0}$ from the constraint (2.8) and consequently $G_{r}>1$. Similarly $a_{0}>0$ when $\gamma<\beta$.

For $\mathbf{Y}_{2}^{*}$ we have $G^{*}=G_{0}$, hence $G_{r}=1$ which implies $A_{11}=0$.

[TABLE]

D.4 Characteristic Equation of the QSP granulopoiesis model

The QSP granulopoiesis model (2.34) has characteristic equation (3.26). From the DDE of the QSP granulopoiesis model (2.34) we compute the linearization matrices and indicate which terms are symmetric and asymmetric. We obtain

[TABLE]

[TABLE]

with

[TABLE]