Processing math: 55%
Research article Special Issues

Equivalences between age structured models and state dependent distributed delay differential equations

  • We use the McKendrick equation with variable ageing rate and randomly distributed maturation time to derive a state dependent distributed delay differential equation. We show that the resulting delay differential equation preserves non-negativity of initial conditions and we characterise local stability of equilibria. By specifying the distribution of maturation age, we recover state dependent discrete, uniform and gamma distributed delay differential equations. We show how to reduce the uniform case to a system of state dependent discrete delay equations and the gamma distributed case to a system of ordinary differential equations. To illustrate the benefits of these reductions, we convert previously published transit compartment models into equivalent distributed delay differential equations.

    Citation: Tyler Cassidy, Morgan Craig, Antony R. Humphries. Equivalences between age structured models and state dependent distributed delay differential equations[J]. Mathematical Biosciences and Engineering, 2019, 16(5): 5419-5450. doi: 10.3934/mbe.2019270

    Related Papers:

    [1] Shaoli Wang, Jianhong Wu, Libin Rong . A note on the global properties of an age-structured viral dynamic model with multiple target cell populations. Mathematical Biosciences and Engineering, 2017, 14(3): 805-820. doi: 10.3934/mbe.2017044
    [2] Cristeta U. Jamilla, Renier G. Mendoza, Victoria May P. Mendoza . Explicit solution of a Lotka-Sharpe-McKendrick system involving neutral delay differential equations using the r-Lambert W function. Mathematical Biosciences and Engineering, 2020, 17(5): 5686-5708. doi: 10.3934/mbe.2020306
    [3] Jordi Ripoll, Jordi Font . Numerical approach to an age-structured Lotka-Volterra model. Mathematical Biosciences and Engineering, 2023, 20(9): 15603-15622. doi: 10.3934/mbe.2023696
    [4] Yijun Lou, Bei Sun . Stage duration distributions and intraspecific competition: a review of continuous stage-structured models. Mathematical Biosciences and Engineering, 2022, 19(8): 7543-7569. doi: 10.3934/mbe.2022355
    [5] Mostafa Adimy, Abdennasser Chekroun, Claudia Pio Ferreira . Global dynamics of a differential-difference system: a case of Kermack-McKendrick SIR model with age-structured protection phase. Mathematical Biosciences and Engineering, 2020, 17(2): 1329-1354. doi: 10.3934/mbe.2020067
    [6] Ming Mei, Yau Shu Wong . Novel stability results for traveling wavefronts in an age-structured reaction-diffusion equation. Mathematical Biosciences and Engineering, 2009, 6(4): 743-752. doi: 10.3934/mbe.2009.6.743
    [7] Qingwen Hu . A model of regulatory dynamics with threshold-type state-dependent delay. Mathematical Biosciences and Engineering, 2018, 15(4): 863-882. doi: 10.3934/mbe.2018039
    [8] Guangrui Li, Ming Mei, Yau Shu Wong . Nonlinear stability of traveling wavefronts in an age-structured reaction-diffusion population model. Mathematical Biosciences and Engineering, 2008, 5(1): 85-100. doi: 10.3934/mbe.2008.5.85
    [9] Fabien Crauste . Global Asymptotic Stability and Hopf Bifurcation for a Blood Cell Production Model. Mathematical Biosciences and Engineering, 2006, 3(2): 325-346. doi: 10.3934/mbe.2006.3.325
    [10] Chang-Yuan Cheng, Shyan-Shiou Chen, Xingfu Zou . On an age structured population model with density-dependent dispersals between two patches. Mathematical Biosciences and Engineering, 2019, 16(5): 4976-4998. doi: 10.3934/mbe.2019251
  • We use the McKendrick equation with variable ageing rate and randomly distributed maturation time to derive a state dependent distributed delay differential equation. We show that the resulting delay differential equation preserves non-negativity of initial conditions and we characterise local stability of equilibria. By specifying the distribution of maturation age, we recover state dependent discrete, uniform and gamma distributed delay differential equations. We show how to reduce the uniform case to a system of state dependent discrete delay equations and the gamma distributed case to a system of ordinary differential equations. To illustrate the benefits of these reductions, we convert previously published transit compartment models into equivalent distributed delay differential equations.


    In Memory of Geoffrey J. Butler and Herbert I. Freedman

    Age structured population models have been used extensively in mathematical biology throughout the past 90 years [1,2] (see [3] for a review). These age structured models describe the progression of individuals through an ageing process by using partial differential equations (PDEs), that can, in certain cases, be reduced to a delay differential equation (DDE) [3,4,5]. When individuals exit the ageing process in a deterministic manner upon reaching a threshold maturation age, the age structured model is typically reduced to a discrete DDE.

    In many populations, the speed at which an individual matures is often only weakly coupled to chronological time and is dynamically controlled by the availability of resources. Consequently, when considering the age of an individual in a population, it is the biological age – and not the chronological age– that is of interest. It is possible to allow for this dynamic accumulation of biological age by including a variable ageing rate in an age structured PDE model. PDE models with variable ageing rates and threshold maturation rates can be reduced to state dependent discrete DDEs. State dependent delays considerably complicate the study of these models, but incorporate external control of the maturation process and increase physiological relevance.

    However, imposing a threshold maturation age does not account for population heterogeneity and implicitly assumes a homogeneous maturation age. Given the importance of individual differences in a population, it is important that intraspecies heterogeneity is included in mathematical models. In light of these observations, we develop a technique to explicitly incorporate maturation age heterogeneity and external control of age accumulation by providing a framework for state dependent distributed DDEs. State dependent distributed DDEs account for a measure of population heterogeneity not present in discrete DDE models while retaining external control of the ageing process. Therefore, distributed DDEs offer a physiologically more realistic manner to model ageing processes in populations [6].

    To derive a state dependent distributed DDE, we consider a general age structured model with a variable ageing rate. We eschew a deterministic maturation process (which would lead to state dependent discrete DDEs), and instead assume that maturation age is a positive random variable A. This random variable defines a density function KA(t) through

    KA(t)=limΔt0P[tAt+Δt]Δt, (1.1)

    which satisfies

    0KA(t)dt=1andKA(t)0t0.

    As shown by Craig et al. [5], Otto and Radons [7], and Bernard [8], replacing existing discrete delays with state dependent delays requires careful attention to how solutions pass across the maturation boundary. Craig et al. [5] derived a "correction" factor to ensure that individuals are not spuriously created or destroyed during maturation. Our work generalises the correction factor derived by Craig et al. [5] for state dependent discrete DDEs to any state dependent DDE. Specifically, our derivation does not rely on a smoothness argument, but arises naturally from the age structured PDE after a careful derivation of the maturation rate.

    We show how the age structured PDE can be reduced to a state dependent distributed DDE. For specific densities KA(t), we show equivalence between the state dependent distributed DDE and state-dependent discrete DDEs with one or two delays or a finite dimensional systems of ordinary differential equations (ODEs). These equivalences arise from the explicit consideration of the ageing process modelled by the distributed DDEs. By applying the linear chain technique to the age variable, instead of the time variable, we are able to establish the desired equivalences. As there is not an available all purpose numerical method capable of solving distributed DDEs, these equivalences allow for the model to be analysed as a DDE and simulated using the highly efficient established techniques for discrete DDEs or ODEs. To illustrate the benefits of the techniques developed here, we consider two previously published models of hematopoietic cell production and show how using distributed DDEs can simplify the analysis of the resulting model.

    The structure of the article is as follows. In Section 2, we study the McKendrick equation for a generic population with a variable ageing rate and random maturation time. By solving the PDE using the method of characteristics, we derive a state-dependent distributed DDE for the general density KA(t) in Theorem 2.1. We discuss the naturally arising correction factor in Section 2.1. To illustrate the benefits of reducing age structured models to DDEs, we show that the resulting DDE preserves non-negativity of initial conditions and perform stability analysis to study the local stability of equilibria in Section 3. By specifying KA(t) to be a degenerate distribution, we recover a state-dependent discrete DDE in Section 4.1. Next, we consider uniform distributions and the equivalent two delay DDE in Section 4.2. In Section 4.3, we study a gamma distributed DDE. Through a generalization of the linear chain technique to include a variable transit rate, we show how this gamma distributed DDE can be reduced to a finite dimensional system of transit compartment ODEs in Section 4.3.1. In Section 5, we formalize the link between variable transit rate compartment models and state dependent delayed processes by converting two previously published transit compartment models to the corresponding distributed DDEs. Finally, we summarize our results with a brief conclusion.

    Consider a population divided into immature and mature compartments in which only mature individuals reproduce. Let n(t,a) denote the number of immature individuals at time t with age a and x(t) denote the number of mature members of the population at time t. The purpose of this section is to establish a state dependent distributed DDE model for x(t).

    We begin with an age structured PDE for the immature population, n(t,a). Immature individuals progress through maturation with a variable ageing rate Va(t), where Va(t) satisfies

    0<VminaVa(t)Vmaxa<.

    Following McKendrick [1], the PDE describing n(t,a) is

    tn(t,a)+Va(t)an(t,a)=[μ(x(t))+h(a)]n(t,a)Va(t)n(t,0)=βx(t)tt0;n(t0,a)=f(a)0a(0,).} (2.1)

    The boundary condition Va(t)n(t,0)=βx(t) that we impose links the creation of immature individuals n(t,0) with the birth rate βx(t). The presence of Va(t) in this boundary term can be understood from the conveyor belt analogy [8,9]. In the following, we assume β>0. The initial conditions n(t0,a)=f(a)0, describes immature individuals with non-zero age at time t0.

    The death rate of immature individuals is given by μ(x(t)) while transition from the immature state to the mature state is modelled by h(a). It is important to note that the transition rate is a function of the age of individuals at time t. Since we expect a link between time and physiological age, we will write a(t). Later, we formalize the weakly coupled relationship between biological and chronological age and justify this notation by finding the characteristics of (2.1).

    We begin by deriving the transition rate from immaturity to maturity, h(a(t)). As mentioned, we assume that the age at which an individual matures is a non-negative random variable A with density function KA(t). The transition rate, h(a(t)), is the instantaneous change in probability that an individual matures at age a(t+Δt), given that the individual has not matured at age a(t). Formally, using the definition of conditional probability,

    h(a(t))=limΔt0P[a(t)Aa(t+Δt) | Aa(t)]Δt=limΔt0P[a(t)Aa(t+Δt)]P[Aa(t)]Δt.

    Multiplying by unity gives

    h(a(t))=1P[Aa(t)]limΔt0P[a(t)Aa(t+Δt)][a(t+Δt)a(t)][a(t+Δt)a(t)]Δt.

    By (1.1) and the derivative of a(t), we obtain

    h(a(t))=KA(a(t))1a(t)0KA(σ)dσddta(t). (2.2)

    The transition (or maturation) rate, h(a(t)), is known as the hazard rate of the random variable A and has applications in modelling failure rates [10,11]. The identical expression for h(a(t)) without considering the conditional maturation probability was derived in [3].

    It is possible that immature individuals create multiple mature individuals upon transitioning to the mature compartment (i.e mitosis), so we model the influx rate into the mature compartment as a function

    F(x(t),0h(s)n(t,s)ds),

    where the integral term

    0h(s)n(t,s)ds (2.3)

    is the number of immature individuals that reach maturity at time t. If mature individuals are cleared at a population dependent rate γ(x(t)), then the mature population satisfies

    ddtx(t)=F(x(t),0h(s)n(t,s)ds)γ(x(t))x(t)x(0)=x0.} (2.4)

    We are now able to establish equivalence between the system of equations describing the populations x(t) and n(t,a) and a distributed DDE. To do this, we partially solve the PDE (2.1) using the method of characteristics.

    Theorem 2.1 (State-Dependent Distributed DDE). Let the immature population n(t,a) satisfy the McKendrick age structured PDE (2.1) with the distribution dependent transition rate h(a(t)) (2.2). Assume that the mature population x(t) is given by (2.4).

    Then, the mature population x(t) satisfies the initial value problem (IVP)

    ddtx(t)=F(x(t),0βx(tφ)Va(t)Va(tφ)exp[ttφμ(x(s))ds]KA(ttφVa(s)ds)dφ)γ(x(t))x(t) (2.5)

    with initial data

    x(s)=ρ(s)s(,t0].

    Proof. The characteristics of equation (2.1) satisfy

    ddφt(φ)=1,andddta(t)=Va(t), (2.6)

    and hence are given by

    t=φ+T0anda(t)=tT0Va(x)dx+a0.

    Along the characteristics, the age structured PDE (2.1) becomes

    ddtn(t,a(t))=[μ(x(t)+KA(a(t))1a(t)0KA(σ)dσVa(t)]n(t,a(t)). (2.7)

    Equation (2.7) has solution

    n(t,a(t))=n(T0,a0)exp[tT0μ(x(s))ds](1a(t)0KA(σ)dσ).

    If a0=0, we use the boundary condition of (2.1) to find

    n(t,a(t))=βx(T0)Va(T0)exp[tT0μ(x(s))ds](1a(t)0KA(σ)dσ), (2.8)

    while, if a0>0, the initial condition of (2.1) gives

    n(t,a(t))=f(a0)exp[tt0μ(x(s))ds](1a(t)0KA(σ)dσ).

    To establish an equivalence between the PDE (2.1) and the distributed DDE (2.5), it is necessary to define suitable initial data x(s)=ρ(s) for s<t0 for the DDE. To do this, it is natural to assume that an an immature individual with positive age a>0 at time t0 was born at sometime s<t0. Since the PDE (2.1) is not defined for s<t0, we are free to prescribe fixed values for Va(s)=Va and μ(x(s))=μ for s<t0. Then, imposing that individuals born at time s<t0 evolved according to the McKendrick Equation, we have a=Va(t0s), or s=t0a/Va. Hence, the initial condition f(a) defines the history function ρ through

    f(a)=βVaρ(t0a/Va)exp[t0t0a/Vaμds]. (2.9)

    Therefore defining x(s)=ρ(s) this way, for s<t0, the solution (2.8) applies.

    Now, we finalize the link between the age structured PDE and the distributed DDE by following the characteristic curves until they intersect with the a=0 axis. Along the characteristic curves, at time t, individuals born at time T0=tφ have age

    at(φ)=tT0Va(x)dx=ttφVa(x)dx

    for φ>0. So we have

    n(t,at(φ))=βx(tφ)Va(tφ)exp[ttφμ(x(s))ds](1at(φ)0KA(σ)dσ).

    At time t, the rate at which individuals mature is

    0h(at(φ))n(t,at(φ))dφ=0KA(at(φ))βx(tφ)Va(t)Va(tφ)exp[ttφμ(x(s))ds]dφ. (2.10)

    By defining, for any density KA(t),

    AK(x(t)):=0KA(at(φ))βx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ, (2.11)

    we have

    0h(at(φ))n(t,at(φ))dφ=Va(t)AK(x(t)).

    Consequently, using (2.11) and defining the history ρ(s) according to (2.9), we have established the equivalence between the system of (2.1) and (2.4) with the distributed DDE (2.5).

    Further inspection of equation (2.10) reveals a ratio of ageing speeds Va(t)/Va(tφ) in the integral term

    0h(at(φ))n(t,at(φ))dφ=0βx(tφ)Va(t)Va(tφ)exp[ttφμ(x(s))ds]KA(σ)dφ.

    The ratio of ageing velocities at the entrance and exit of the ageing process acts as a correction factor. As shown by Bernard [8] and Craig et al. [5], models without the correction factor allow for spurious creation of individuals during maturation and some state-dependent DDEs have missed this important correction factor. Solutions of models without this correction factor do not necessarily preserve nonnegativity of initial data [8].

    Craig et al. [5] derived the correction factor by carefully accounting for the number of cells crossing the maturation threshold in a discrete state-dependent DDE. In discrete DDEs, individuals mature following a deterministic process after accruing a specific threshold age, so the maturation boundary is well-defined. The derivation of the correction factor was based on the smoothness of the solution crossing the fixed maturation boundary. However, the idea of a fixed maturation boundary does not extend to random maturation ages. Consequently, the derivation of the correction factor by Craig et al. [5] does not generalise to generic distributed DDEs.

    Our derivation of the state-dependent distributed DDE produces the same correction factor through the instantaneous maturation probability, h(a(t)). The derivation of h(a(t)) in equation (2.2) produces the term Va(t) by accounting for the change of maturation probability due to the variable accumulation of age at time t. For a degenerate distribution, as shown in Section 4.1, we obtain precisely the same ratio as [5].

    Replacing an age structured PDE by a DDE eliminates the need to explicitly model the ageing populations, which can be difficult to measure experimentally. DDEs offer a natural framework that explicitly incorporates delays and identifies the relationship between the current and past states. This can facilitate communication between mathematical biologists and biologists and physiologists. In particular, the explicit presence of the delay term allows for simple calculation of mean delay time. As shown in [12], models of delayed processes without DDEs do not always accurately calculate the mean delay time. However, DDEs typically define infinite dimensional semi-dynamical systems, which can introduce mathematical difficulties.

    As we have seen in Theorem 2.1, partially solving an age structured PDE may lead to a DDE. As such, analysing these partially solved systems can be simpler than studying the corresponding PDE. As an example, we analyse the state-dependent distributed DDE in equation (2.5). Define

    ˉx(t)=Va(t)AK(t)=Va(t)0KA(at(φ))βx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ, (3.1)

    and consider the IVP

    ddtx(t)=F[x(t),ˉx(t)]γ(x(t))x(t)t>t0x(s)=ρ(s)s(,t0],} (3.2)

    where F(x,y)C1(R2,R) and γ(x(t))C1(R,R) with

    F(x,y)>0ifx>0ory>0,F(0,0)=0,andγ(x(t))<γmax<. (3.3)

    We recall that A is the random variable representing the maturation age of immature individuals. The history function, ρ(s), is chosen to belong to the space L1(A) where

    KA(t)=dAdλ,

    and λ is the Lebesgue measure on R. L1(A) satisfies the axioms for a phase space given in [13] and [14], so the solution of the IVP (3.2) exists and is unique in L1(A). In population modelling, it is likely that any realistic history is uniformly continuous and bounded. The space of bounded and uniformly continuous functions is a subspace of L1(A) and is a suitable phase space.

    The age structured PDE (2.1) describes population dynamics in the presence of a maturation time. Consequently, solutions of (3.2) must represent a population, and in particular, remain non-negative. However, the presence of delays in other models may lead to solutions that do not remain non-negative, as noted in [15]. We begin our analysis by showing that the solution of the IVP (3.2), x(t), evolving from non-negative initial conditions remains non-negative. This property is a natural requirement for models of population dynamics.

    Proposition 3.1. Let F(x,y) and γ(x(t)) satisfy equations (3.3). Moreover, assume that the history function satisfies

    ρ(s)0s(,t0].

    Then, the solution of the IVP (3.2) remains non-negative for all time t>t0.

    Proof. As ρ(s)0, it is simple to see that

    ˉx(t0)=Va(t0)0KA(at0(φ))βρ[t0φ]Va(t0φ)exp[t0t0φμ(x(s))ds]dφ0.

    We have a series of cases.

    1) If ρ(t0)=x(t0)>0 then F(x(t0),ˉx(t0))>0. Therefore,

    ddtx(t)=F(x(t),ˉx(t))γ(x(t))x(t)γ(x(t))x(t)>γmaxx(t)

    and using Gronwall's inequality, we have

    x(t)ρ(t0)exp(γmax[tt0])>0.

    2) If ρ(t0)=0 and ρ(s)=0 A-almost everywhere in (,t0), then x(t)=0 is the solution of the IVP.

    3) Finally, if ρ(t0)=0 and ρ(s)>0 on a set of A-positive measure in (,t0) then ˉx(t0)>0 and

    ddtx(t)|t=t0=F(x(t0),ˉx(t0))γ(t0)x(t0)=F(0,ˉx(t0))>0.

    Consequently, x(t) becomes positive immediately and Case 3 reduces to Case 1.

    Therefore, solutions of the IVP (3.2) remain non-negative for all time t>t0.

    We continue the analysis of equation (2.5) by studying the local stability of equilibrium solutions. To do this, let x(t)=xL1(A) be an equilibrium of the IVP (3.2), so

    F(x,ˉx)=γ(x)x. (3.4)

    At the equilibrium x(t)=x, Va(t)=Va, so

    at(φ)=ttφVa(s)ds=VaφandVa(t)Va(tφ)=1

    then, by evaluating (3.1) with x(t)=x, the homeostatic delayed term ˉx in (3.4) satisfies

    ˉx=0βxKA(Vaφ)exp[μφ]dφ=βVaxL[KA](μ/Va),

    where L[f](s) is the Laplace transform of f(x) evaluated at s.

    Hence, ˉx is a function of the density KA(t). However, if desired, it is possible to vary the homeostatic death rate μ to ensure that the equilibria value x does not change for different densities KA(t) as shown in [6].

    Set z(t)=x(t)x, and for z(t) small– similar to the discrete state dependent delay case considered in [16]– freeze the ageing and clearance rates at their homeostatic values, so Va(t)=Va and μ(s)=μ. By expanding exponential integral to leading order following [6], it is possible to relax the assumption that μ(s)=μ. Now, define ˉz(t)=ˉx(t)ˉx so that

    ˉz(t)=0KA(Vaφ)βx[tφ]exp[μφ]βxKA(Vaφ)exp[μφ]dφ=0KA(Vaφ)βz[tφ]exp[μφ]dφ, (3.5)

    and the equilibrium is translated to the origin. Then, the differential equation for z(t) is

    ddtz(t)=F(z(t)+x,ˉz(t)+ˉx)γ(x(t))z(t)γ(x(t))x.

    By making the ansatz

    z(t)=Ceλt,

    we compute the expression for ˉz(t) from (3.5)

    ˉz(t)=Cz(t)0KA(Vaφ)βeλφ[exp[μφ]]dφ=Cz(t)βVaL[KA]([μ+λ]/Va).

    Therefore

    ddtz(t)=k1z(t)+k2βVaL[KA]([μ+λ]/Va)z(t)γz(t)+O(z2)

    where γ=xγ(x(t))|x=x, k1=aF(a,b)|(x,ˉx) and k2=bF(a,b)|(x,ˉx). Dropping nonlinear terms, the linearised equation is

    ddtz(t)=(k1γ)z(t)+k2βVaL[KA]([μ+λ]/Va)z(t). (3.6)

    The characteristic equation corresponding to (3.6) is

    0=λ(k1γ)k2βVaL[KA]([μ+λ]/Va). (3.7)

    Through a standard analysis, we study the local stability of the equilibrium x for a density KA(t).

    Proposition 3.2. 1) If

    |k2|βVaL[KA](μ/Va)<γk1,

    the equilibrium point x is locally asymptotically stable.

    2) If

    k2βVaL[KA](μ/Va)>γk1,

    the equilibrium point x is unstable.

    Proof. 1) Let λ be a root of (3.7) and assume for contradiction that (λ)0. We necessarily have

    λ=(k1γ)+k2βVaL[KA]([μ+λ]/Va),

    and we calculate

    (λ)=(k1γ)+k2βVa[L[KA]([μ+λ]/Va)].

    We note that

    k2βVa[L[KA]([μ+λ]/Va)]|k2βVaL[KA]([λ+μ]/Va)|.

    While, for arbitrary ν=νr+iνiC,

    |k2βVaL[KA]([μ+ν]/Va)|=|k2|βVa|0exp[(μ+νr+iνi)φ]KA(Vaφ)dφ||k2|βVa0exp[(μ+νr)φ]KA(Vaφ)|eiνiφ|dφ=|k2|βVaL[KA]([μ+νr]/Va).

    Moreover, if νr0,

    |k2|βVaL[KA]([μ+νr]/Va)|k2|βVaL[KA](μ/Va).

    Therefore, using the assumption in 1), we find

    (λ)=(k1γ)+k2βVa[L[KA]([λ+μ]/Va)](k1γ)+|k2|βVaL[KA](μ/Va)<0,

    which is a contradiction, so no such λ can exist. Therefore, all roots of the characteristic equation have negative real part and the equilibrium is stable.

    2) To show instability, we will prove that there must be one characteristic root with positive real part. Define

    g(λ):=k1γλ+k2βVaL[KA]([λ+μ]/Va),

    and note that g(λ) is continuous with

    g(0)=k1γ+k2βVaL[KA](μ/Va)>0andlimλg(λ)=.

    Then, there must be a real λ>0 such that g(λ)=0. The equilibrium is therefore unstable.

    We note that if k2>0, i.e. the production of mature individuals is controlled through positive feedback with the number of maturing individuals at time t, then Proposition 3.2 completely characterizes the local stability of x. If k2<0, it seems likely that x would lose stability through a Hopf bifurcation, similar to the discrete delay case. A similar analysis was done in the constant ageing rate in [17]. However, Yuan and Bélair [17] did not consider death of immature individuals, nor the linear clearance of mature individuals which corresponds to μ=γ=0.

    Next, we study the DDE found in Theorem 2.1 for various density functions. By first considering the characteristic equation (3.7) for specific densities KA(t), we motivate the reduction of these population models to familiar discrete DDEs and transit compartment ODEs. In the discussion that follows, we once again assume that xL1(A) is an equilibrium point so that μ(x)=μ and Va(t)=Va. Denote the homeostatic maturation time as the first moment of the random variable A with constant ageing rate Va,

    τ=0tKA(Vat)dt.

    Consequently, the expected homeostatic maturation age is given by T=Vaτ.

    We first consider the degenerate distribution concentrated at T and recover the familiar state dependent discrete DDE. Next, we use a linear chain-type technique to reduce state dependent uniformly distributed DDEs to a system involving two state dependent delays. Finally, we show how to reduce a gamma distributed DDE to a transit compartment system of ODEs.

    However, true equivalence between the distributed DDE and the reduced form does not follow directly. We must take care when prescribing initial conditions and history functions so that solutions of the different formulations are in fact equivalent. Only then do these reductions allow for the use of the highly efficient numerical methods available for discrete DDEs and ODEs available in most programming languages.

    Assuming that maturation is a deterministic process and occurs after achieving the threshold age T implies that KA(t) is the degenerate distribution concentrated at T with

    KA(ttφVa(s)ds)=δ(ttφVa(s)dsT). (4.1)

    where δ(x) is the Dirac delta function. In the deterministic case, all individuals mature at precisely the same age T. At the equilibrium x, using (3.7), the characteristic equation is

    0=λ(k1γ)k2βVaexp[(μ+λ)T/Va]=λ(k1γ)k2βVaexp[(μ+λ)τ], (4.2)

    which is exactly the characteristic equation of a discrete DDE. This is unsurprising, since it is well known that threshold conditions lead to discrete DDEs [4,7].

    Returning to the DDE (2.5) with KA(t) given by (4.1), the threshold maturation age T allows us to calculate when an individual that matures at time t began maturation. The maturation time, τ(x(t)), must satisfy the implicit threshold condition

    T=ttτ(x(t))Va(s)ds. (4.3)

    We use the definition of τ(x(t)) to evaluate the convolution integral given in (2.11) to find

    Aδ(t)=0δ(ttφVa(s)dsT)βx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ=βx[tτ(x(t))]Va(tτ(x(t)))exp[ttτ(x(t))μ(x(s))ds].

    Consequently, the corresponding IVP to (2.5) with state dependent discrete delay is

    ddtx(t)=F(x(t),βx[tτ(x(t))]exp[ttτ(t)μ(x(s))ds]Va(t)Va(tτ(t)))γx(t)x(s)=ρ(s)s(,t0].} (4.4)

    To implement (4.4) numerically, it is necessary to solve (4.3) to find the maturation time τ(x(t)). This can be done by differentiating (4.3) to find

    ddtτ(x(t))=1Va(t)Va(tτ(x(t))), (4.5)

    and imposing the correct initial condition so that the solution of (4.5) also solves (4.3).

    In the case that ρ(s)=x, then it is simple to set τ(0)=τ. However, for more general initial data ρ(s), choosing an appropriate initial condition for (4.5) can be delicate [7].

    Then, we can solve the discrete state dependent DDE by solving the system of equations given by (4.4) and (4.5). Hence the age structured PDE framework in Section 2 offers an alternative to the "moving threshold" method to derive state dependent DDEs as described in [7].

    We consider uniformly distributed DDEs centered about the expected homeostatic maturation age T. In the simplest case, the uniform distribution defines lower and upper threshold ages and assigns equal weight to each age falling between the thresholds. The probability density function corresponding to a uniform distribution centred at T is

    KU(a)={12Vaδifa[TVaδ,T+Vaδ]0otherwise. (4.6)

    At the equilibrium x, with KA(t) given by the uniform density (4.6), the characteristic equation (3.7) is

    0=λ(k1γ)k2(βVa)(12δVa[λ+μ]/Va)[e(λ+μ)(TVaδ)/Vae(λ+μ)(T+Vaδ)/Va]=λ(k1γ)k2(βVa)(12δ(λ+μ))[e(λ+μ)(τδ)e(λ+μ)(τ+δ)]. (4.7)

    TVaδ and T+Vaδ represent the minimal and the maximal ages at which an individual can mature. Due to the variable ageing rate, the minimal and maximal delay times, τmin(x(t)) and τmax(x(t)), are state dependent, and implicitly defined by

    TVaδ=ttτmin(x(t))Va(s)dsandT+Vaδ=ttτmax(x(t))Va(s)ds.

    We note that, at homeostasis, Va(s)=Va so

    TVaδ=τmin(x)VaandT+Vaδ=τmax(x)Va.

    Recalling that T=Vaτ, the terms τδ and τ+δ in (4.7) correspond to the minimal and maximal homeostatic delay times.

    The presence of minimal and maximal delay terms in (4.7) hints that a uniformly distributed DDE may be reducible to a discrete DDE with two distinct delays.

    Inserting the uniform density (4.6) into the convolution integral (2.11) gives

    AU(t)=0KU(ttφVa(s)ds)βx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ=τmax(t)τmin(t)12Vaδβx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ.

    Thus the state dependent uniform distributed DDE is

    ddtx(t)=F(x(t),AU(t)Va(t))γ(x(t))x(t)x(s)=ρ(s),s(,t0].} (4.8)

    Next, we show that (4.8) can be reduced to an IVP with two state dependent discrete delays. Once again, this is advantageous, as numerical algorithms for systems of state dependent discrete DDEs are available in most programming languages.

    We begin by formalizing the link between uniformly distributed DDEs and discrete DDEs that was hinted at in (4.7). To do this, we show to write the delay kernel as the solution of a differential equation. This strategy is also used in the well known linear chain technique (see [18]), which we generalise in Section 4.3.1. However, unlike the linear chain technique, we will not recover a system of ODEs, but rather a system of differential equations with two state dependent discrete delays. The technique here can also be adapted to "tent" like distributions (see [19]).

    Lemma 4.1. AU(t) satisfies the differential equation

    ddtAU(t)=12Vaδ[βx[tτmin(t)]Va(tτmin(t))exp[ttτmin(t)μ(x(s))ds]Va(t)Va(tτmin(t))βx[tτmax(t)]Va(tτmax(t))exp[ttτmax(t)μ(x(s))ds]Va(t)Va(tτmax(t))]μ(x(t))AU(t). (4.9)

    Proof. Similar to the linear chain technique, we differentiate AU(t) using Leibniz's rule to find

    ddtAU(t)=12Vaδ[βx[tτmax(t)]Va(tτmax(t))exp[ttτmax(t)μ(x(s))ds]ddtτmax(t)βx[tτmin(t)]Va(tτmin(t))exp[ttτmin(t)μ(x(s))ds]ddtτmin(t)]+12δτmax(t)τmin(t)ddt(βx(tφ)Va(tφ)exp[ttφμ(x(s))ds])dφ.

    We note that

    ddt(βx(tφ)Va(tφ)exp[ttφμ(x(s))ds]dφ)=ddφ(βx(tφ)Va(tφ)exp[ttφμ(x(s))ds])μ(x(t))βx(tφ)Va(tφ)exp[ttφμ(x(s))ds],

    so that, integrating by parts,

    τmax(t)τmin(t)ddt(βx(tφ)Va(tφ)exp[ttφμ(x(s))ds])dφ=(12δβx(tφ)Va(tφ)exp[ttφμ(x(s))ds])|τmax(t)φ=τmin(t)μ(x(t))AU(t).

    Consequently, the derivative of AU(t) is

    ddtAU(t)=12Vaδ[βx[tτmax(t)]Va(tτmax(t))exp[ttτmax(t)μ(x(s))ds](ddtτmax(t)1)βx[tτmin(t)]Va(tτmin(t))exp[ttτmin(t)μ(x(s))ds](ddtτmin(t)1)]μ(x(t))AU(t).

    To finish the proof, we note that, similar to (4.5), τmin(x(t)) and τmax(x(t)) solve the following differential equations

    ddtτmin(x(t))1=Va(t)Va(tτmin(x(t)))andddtτmax(x(t))1=Va(t)Va(tτmax(x(t))). (4.10)

    The identities in equation (4.10) give (4.9).

    By writing the delay term AU(t) as a solution of a differential equation, we are able to reduce the distributed DDE to a system with state dependent discrete delays. Once again, this allows for simulation of the distributed DDE (4.8) using existing techniques. This relationship is formalized in the following theorem.

    Theorem 4.2. The IVP (4.8) is equivalent to the IVP with the following system of discrete delay differential equations

    ddtx(t)=F(x(t),y(t)Va(t))γ(x(t))x(t)ddty(t)=12δ[βx[tτmin(t)]ˆVa(tτmin(t))exp[ttτmin(t)μ(x(s))ds]Va(t)Va(tτmin(t))βx[tτmax(t)]ˆVa(tτmax(t))exp[ttτmax(t)μ(x(s))ds]Va(t)Va(tτmax(t))]μ(x(t))y(t).} (4.11)

    with suitably chosen initial data.

    Proof. Using Lemma 4.1, it is simple to see that

    \begin{equation} y(t)V_a(t) = A_U(t)V_a(t), \end{equation} (4.12)

    and the other terms in the differential equations are identical if the initial data are equivalent. It therefore remains to show that we can choose suitable history functions for the distributed and discrete DDEs. For the history function of the distributed DDE (4.8), \rho(s) , setting the initial data of (4.11) to be

    \begin{equation*} x(s) = \rho(s) \end{equation*}

    and

    \begin{equation*} y(t_0) = \int_{\tau_{min}(t)}^{\tau_{max}(t)} \frac{1}{2\delta} \frac{ \beta \rho(t_0- \varphi) }{\hat{V}_a (t_0- \varphi)} \exp \left[ - \int_{t_0- \varphi}^{t_0} \mu(\rho(s)) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation*}

    gives the desired equivalence [19]. To convert from (4.11) with history function x(s) = \eta(s) to (4.8), y(t_0) must satisfy

    \begin{equation} y(t_0) = \int_{\tau_{min}(t)}^{\tau_{max}(t)} \frac{1}{2\delta} \frac{ \beta \eta(t_0- \varphi) }{\hat{V}_a (t_0- \varphi)} \exp \left[ - \int_{t_0- \varphi}^{t_0} \mu(\eta(s)) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation} (4.13)

    By taking the initial data for (4.8) to be x(s) = \eta(s) , we see that this condition is sufficient for equivalence of (4.11) and (4.8). Now, if (4.13) does not hold, then (4.12) cannot be satisfied at t = t_0 , so (4.13) is a necessary and sufficient condition to be able to convert the system of DDEs (4.11) into the distributed DDE (4.8) with x(s) = \eta(s) .

    Finally, we study gamma distributed DDEs and show how to reduce the state-dependent gamma distributed DDE to a transit chain of ODEs. The probability density function of the gamma distribution is

    \begin{equation} g_b^j(x) = \frac{b^jx^{j-1}e^{-bx}}{\Gamma(j)}, \end{equation} (4.14)

    where j, b \in \mathbb{R} .

    Again, let \mathcal{T} denote the mean maturation age and fix j > 0 . Then, we have the following relationships

    \begin{equation*} \mathcal{T} = j/V_a^*, \quad \sigma^2 = j/(V_a^*)^2, \quad {\rm and} \quad K_{g}(\sigma) = g_{V_a^*}^j(\sigma), \end{equation*}

    where \sigma^2 is the variance of the gamma distribution and we set b = V_a^* .

    Calculating (3.7) for the gamma density in (4.14) gives

    \begin{equation} 0 = k_1-\gamma - \lambda + k_2 \left( \frac{\beta}{V_a^*} \right) \left( \frac{(V_a^*)^j}{(V_a^*+[\lambda+\mu^*]/V_a^*)^j} \right) . \end{equation} (4.15)

    Now, we use the relationships V_a^* = j/\mathcal{T} and \mathcal{T} = \tau^* V_a^* to rewrite the characteristic function as

    \begin{align*} k_1-\gamma - \lambda + k_2 \left( \frac{\beta}{V_a^*} \right) \left( \frac{1}{(1+\frac{\lambda+\mu^*}{V_a^*}^2)^j} \right)& = k_1-\gamma - \lambda + k_2 \left( \frac{\beta}{V_a^*} \right) \left( \frac{1}{(1+\frac{\mathcal{T}(\lambda+\mu^*)}{ V_a^* j})^j} \right) \\ & = k_1-\gamma - \lambda + k_2 \left( \frac{\beta}{V_a^*} \right) \left( \frac{1}{(1+\frac{\tau^*(\lambda+\mu^*)}{j})^j} \right). \end{align*}

    Using a common denominator gives

    \begin{equation} 0 = \left(k_1-\gamma-\lambda\right)\left(1+\frac{\tau^*(\lambda+\mu^*)}{j}\right)^j+k_2\frac{\beta}{V_a^*}. \end{equation} (4.16)

    Now, we consider multiple cases for the parameter j . If j \in \mathbb{N} , then (4.16) is a polynominal of degree j+1 , with j+1 roots. This is markedly different than the generic distributed DDE, as the characteristic equation (3.7) is typically a transcendental function of \lambda with infinitely many characteristic values. Now, with j = n/m \in \mathbb{Q} , we can rearrange (4.16) to

    \begin{equation*} \left(k_1-\gamma-\lambda\right)\left(1+\frac{\tau^*(\lambda+\mu^*)}{j}\right)^j = -k_2\frac{\beta}{V_a^*} , \end{equation*}

    and raising both sides of the equality to the power m gives

    \begin{equation} 0 = (k_1-\gamma-\lambda)^m(1+\frac{\tau^*(\lambda+\mu^*)}{j})^n + \left(-k_2\frac{\beta}{V_a^*} \right)^m . \end{equation} (4.17)

    Not all solutions of (4.17) will necessarily satisfy (4.15). However, every solution of (4.15) will satisfy (4.17). Moreover, (4.17) is a polynomial with m+n roots, so (4.15) with j = n/m \in \mathbb{Q} has at most m+n roots. However, if the parameter j is not rational, then (4.16) is once again a transcendental equation with possibly infinitely many roots.

    The relationship between the number of characteristic values and the parameter j leads to interesting questions. If j \in \mathbb{N} increases by unit steps, then the characteristic equation gains precisely one root. However, if j increases smoothly between j and j+1 , do characteristic values spring in and out of existence depending on the rationality of j ? This question, while important, is outside the scope of the current work.

    Having studied the characteristic equation of gamma distributed DDEs, we proceed to write down the gamma distributed DDE. We have parametrized the gamma distribution so that at homeostasis, the mean delay time is \tau^* . The variable ageing velocity must then be scaled so that at homeostasis, individuals age chronologically. Therefore, we define the scaled ageing velocity

    \begin{equation} \hat{V}_a(t) = \frac{V_a(t)}{V_a^*}, \end{equation} (4.18)

    and will use \hat{V}_a(t) throughout the remainder of our study. The scaled density function g_{V_a^*}^j(a_t(\varphi)) is given by

    \begin{equation*} g_{V_a^*}^j\left( \int_{t- \varphi}^t \hat{V}_a(s) \mathrm{d} s\right) = \frac{(V_a^*)^j}{\Gamma(j)}\left[ \int_{t- \varphi}^t \hat{V}_a(s) \mathrm{d} s\right]^{j-1} \exp \left[-V_a^* \int_{t- \varphi}^t \hat{V}_a(s) \mathrm{d} s \right]. \end{equation*}

    By inserting g_{V_a^*}^j(a_t(\varphi)) into equation (2.11), we define

    \begin{equation} A_{g}(t) = \int_{0}^{\infty} g_{V_a^*}^j \left( \int_{t- \varphi}^t \hat{V}_a(s) \mathrm{d} s\right) \frac{ \beta x(t- \varphi) }{\hat{V}_a (t- \varphi)} \exp \left[ - \int_{t- \varphi}^t \mu(x(s)) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation} (4.19)

    Then, the IVP with a state-dependent distributed DDE corresponding to equation (2.5) is

    \begin{equation} \left. \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} x(t) & = F \left[ x(t), V_a(t) A_{g}(t)\right] - \gamma(x(t)) x(t)\\ x(s) & = \rho(s), s \in (-\infty, t_0]. \end{aligned} \right \} \end{equation} (4.20)

    As we show in Section 5, equivalent models to (4.20) have been used in pharmacokinetic modelling. However, these models typically take the form of finite dimensional systems of ODEs and the direct link between these ODEs with variable transit rates and (4.20) has not been established previously.

    The finitely many roots of equation (4.15) for integer j\in \mathbb{N} suggest that there is a finite dimensional representation of the DDE (4.20). The link between gamma distributed DDEs and transit chain ODEs with constant transit rates has been known since at least 1961 [20]. The method entered into the English literature in the works of MacDonald [21] as the linear chain trick or the linear chain technique.

    Just as in Section 4.2, the linear chain technique consists of replacing the convolution integral (4.19) by the solution of a system of differential equations. To do this, we will exploit the fact that, for j\in \mathbb{N} ,

    \begin{equation} \frac{ \mathrm{d}}{ \mathrm{d} x} g_b^1(x) = -b g_b^1(x) \quad {\rm and} \quad \frac{ \mathrm{d}}{ \mathrm{d} x} g_b^j(x) = b[g_b^{j-1}(x)-g_b^j(x)]. \end{equation} (4.21)

    The linear chain technique has been used extensively in pharmacology to model delayed drug absorption and action. However, typical applications of the technique require that transition rates between compartments are constant and identical. De Souza et al. [12] developed an adapted linear chain technique that allows for variable transition rates by rescaling time in a non-linear way. This non-linear time rescaling leads to difficulties in establishing a link between time rescaled simulations and time series patient data [12]. Here, we provide an alternative technique that allows for variable transition rates between compartments without rescaling time.

    We first show how to write (4.19) as the solution of a system of ordinary differential equations.

    Lemma 4.3. For j\in \mathbb{N} , A_g(t) = x_j(t) where \{ x_i(t) \}_{i = 1}^j satisfies

    \begin{equation} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} x_1(t) & = \frac{ \beta x(t) }{\hat{V}_a (t)} - V_a(t)x_1(t) -\mu(x(t)) x_1(t)\\ \notag \frac{ {\rm d}}{ {\rm dt}} x_i(t) & = V_a(t)\left[ x_{i-1}(t) - x_{i}(t)\right] - \mu(x(t)) x_i(t)\quad {\rm for} \quad i = 2, 3, ..., j. \notag \end{aligned} \label{Eq:TransitChainODEGeneral} \end{equation}

    Proof. We first note that

    \begin{equation*} g_{V_a^*}^i\left( \int_{t}^t \hat{V}_a(s) \mathrm{d} s \right) = \left \{ \begin{array}{lll} V_a^* & {\rm if} & i = 1 \\ 0 & {\rm if} & i = 2, 3, ..., j. \\ \end{array} \right. \end{equation*}

    Then using (4.21) and (2.6), the chain and Leibniz rules show that

    \begin{equation*} \frac{ {\rm d}}{ {\rm dt}} g_{V_a^*}^1 \left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s\right) = -V_a^*\hat{V}_a(t) g_{V_a^*}^1 \left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s\right) = -V_a(t) g_{V_a^*}^1\left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s\right) \end{equation*}

    while, for i = 2, 3, 4, ... ,

    \begin{equation*} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} g_{V_a^*}^i \left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s \right) & = V_a^* \hat{V}_a(t) \left[ g_{V_a^*}^{i-1}\left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s \right) -g_{V_a^*}^i\left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s \right) \right]\\ & = V_a(t) \left[ g_{V_a^*}^{i-1} \left( \int_{ \varphi}^t \hat{V}_a(s) \mathrm{d} s\right)-g_{V_a^*}^i\left( \int_{ \varphi }^t \hat{V}_a(s) \mathrm{d} s\right) \right]. \end{aligned} \end{equation*}

    Now we define,

    \begin{equation*} a_t(x) = \int_{t-x}^t \hat{V}_a(s) \mathrm{d} s \end{equation*}

    and, for i = 1, 2, ..., j ,

    \begin{equation} x_i(t) = \int_{-\infty}^{t} g_{V_a^*}^i (a_t(t- \varphi)) \frac{ \beta x( \varphi) }{V_a ( \varphi)} \exp \left[ - \int_{ \varphi}^t \mu(x(s)) \mathrm{d} s\right] \mathrm{d} \varphi, \end{equation} (4.22)

    and note that, after making the change of variable u = t- \varphi in A_g(t) ,

    \begin{equation*} x_j(t) = \int_{-\infty}^{t} g_{V_a^*}^j (a_t(t- \varphi)) \frac{ \beta x[ \varphi] }{V_a ( \varphi)} \exp \left[ - \int_{ \varphi}^t \mu(x(s)) \mathrm{d} s\right] \mathrm{d} \varphi = A_g(t). \end{equation*}

    Now, by differentiating (4.22) using the Leibniz rule, the transit chain x_i(t) satisfies the system of equations

    \begin{align} \frac{ {\rm d}}{ {\rm dt}} x_1(t) & = \frac{ \beta x(t) }{\hat{V}_a (t)} - V_a(t)x_1(t) -\mu(x(t)) x_1(t)\\ \frac{ {\rm d}}{ {\rm dt}} x_i(t) & = V_a(t)\left[ x_{i-1}(t) - x_{i}(t)\right] - \mu(x(t)) x_i(t)\quad {\rm for} \quad i = 2, 3, ..., j. \end{align} (4.23)

    Importantly, Lemma 4.3 ensures that

    \begin{equation} V_a(t) A_{g}(t) = V_a(t) x_j(t). \end{equation} (4.24)

    Now, we can use the relationship between equations (4.22) and (4.24) to establish the following theorem:

    Theorem 4.4 (Finite Dimensional Representation). The distributed state dependent DDE (4.20) with j \in \mathbb{N} is equivalent to the finite dimensional transit compartment ODE system given by

    \begin{equation} \left. \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} x(t) & = F(x(t), V_a(t)x_j(t) )- \gamma(x(t)) x(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} x_1(t) & = \frac{ \beta x(t) }{\hat{V}_a (t)} - V_a(t)x_1(t)-\mu(x(t)) x_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} x_i(t) & = V_a(t)\left[ x_{i-1}(t) - x_{i}(t)\right] -\mu(x(t)) x_i(t) \quad {\rm for} \quad i = 2, 3, ..., j. \end{aligned} \right \} \end{equation} (4.25)

    Proof. Lemma 4.3 ensures that the differential equations are equivalent. Therefore, we need only construct appropriate initial data for the distributed DDE and ODE formulation. For a history function \rho(s) of (4.20), we set, for i = 1, 2, ..., j,

    \begin{equation} x_i(0) = \int_{-\infty}^{0} g_{V_a^*}^i (a(- \varphi)) \frac{ \beta \rho( \varphi) }{V_a ( \varphi)} \exp \left[ - \int_{ \varphi}^t \mu(\rho(s)) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation} (4.26)

    If \mu(s) = \mu^* is constant and the initial conditions satisfy

    \begin{equation*} x_i(0) = \left(\frac{V_a^*}{V_a^*+\mu^*}\right)^i x_1(0), \end{equation*}

    it is simple to choose \rho(s) = x_1(0) . However, in the more general case with \mu(t) \neq \mu^* and arbitrary ODE initial conditions x_i(0) = \alpha_i of (4.22), we can use a similar method to Cassidy and Humphries [6] to construct one of the infinitely many appropriate history functions.

    A form of the expression for the variable age transit chain in equation (4.22) was derived in [22] to study the equivalence between lifespan and transit compartment models in pharmacodynamics. However, the derivation did not include the underlying age structured PDE and was specific to the gamma distribution. Gurney et al. [23] derived a similar expression for the density of individuals progressing through a specific stage of maturation from a balance equation. However, they did not explicitly formulate the underlying DDE nor did they derive the correct initial conditions for each of the transit compartments. Consequently, they did not show equivalence between the transit compartment formulation and the DDE.

    Remark 4.5 (Recipe for equivalency between ODEs and gamma distributed DDEs) We note that the finite dimensional representation of (4.20) with j\in \mathbb{N} includes a transit compartment chain. Due to the equivalence between (4.20) and (4.24), we are able to identify the ingredients needed to transform a transit compartment ODE such as (4.24) into a DDE such as (4.20). We first consider

    \begin{equation*} \frac{ {\rm d}}{ {\rm dt}} x_1(t) = \frac{ \beta x(t) }{\hat{V}_a (t)} - V_a(t)x_1(t)-\mu(x(t)) x_1(t). \end{equation*}

    From the equation for x_1(t) , we can easily identify the ratio \beta x(t)/\hat{V}_a (t) as the rate at which individuals in the 1st compartment are created. Next, by considering the rate at which individuals enter the second compartment,

    \begin{equation} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} x_1(t) & = \frac{ \beta x(t) }{\hat{V}_a (t)} - V_a(t)x_1(t)-\mu(x(t)) x_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} x_2(t) & = V_a(t)x_1(t) - V_a(t)x_{2}(t) -\mu(x(t)) x_2(t), \end{aligned} \end{equation} (4.27)

    we find the (possibly variable) transit rate between compartments. Then, a process of elimination immediately yields the mortality rate \mu(x(t)) (if \mu(x(t)) < 0 , then population growth rather than decay is occurring through the transit chain). The creation and transit rates also yield the homeostatic ageing rate via (4.18). Further inspection of (4.19) shows that these rates are all that are needed to transform the transit chain ODE to a distributed DDE.

    We note that the classic linear chain technique (see [18]) is a special case of Remark 4.5 where the ageing velocity, V_a(t) , is constant.

    Sometimes, analysis of distributed DDEs is more tractable and simpler than that of a high dimensional equivalent ODE system. For example, by rescaling time, de Souza et al. [12] converted Quartino's ODE transit compartment model of granulopoiesis into a distributed DDE [24]. The distributed DDE formulation proved to be much more analytically tractable than the ODE case, and was used to show the positivity of solutions and establish the local stability of equilibrium solutions.

    However, due to the lack of a general numerical algorithm, simulation of distributed DDEs must be handled on a case by case basis. Simulation of transit compartment ODEs is routine in many programming languages and can be used for the calibration of models to existing data. Once calibrated, mathematical models can be simulated and used in a predictive manner. Consequently, by converting models between the equivalent distributed DDE or ODE formulations, researchers can use the form of the model that is most suitable to their needs.

    The hematopoietic system controls blood cell production and, through tight cytokine control, is able to quickly respond to challenges, including infection and blood loss. Cytokines control hematopoietic output by varying effective proliferation and maturation rates in each hematopoietic lineage. As cells are not produced instantaneously, there is necessarily a delay between cytokine signal and production response. Mathematical models have been used to understand the complex dynamics observed in so-called dynamical diseases since the 1970s [25,26,27]. Existing mathematical models of hematopoiesis have included discrete, distributed and state-dependent DDEs [5,9,28,29,30] as well as transit compartment models [31,32,33].

    Here, we use the equivalence between state dependent distributed DDEs and ODE transit compartment models derived in Section 4.3.1 to convert two previously published ODE models of hematopoietic cell production to their equivalent state-dependent distributed DDEs. The ODE models specify the entrance rate of individuals into the maturation compartment and the maturation speed, V_a(t) , which allows for the calculation the birth rate \beta of immature individuals. As these models involve more than one population, the birth rate \beta is no longer constant but is a function of other populations in the model.

    In the first example, we show how a model of reticulocyte production can be reduced to a renewal equation whose dynamics are completely characterized by a simple system of ordinary differential equations.

    In the second example, we extend the framework of Section 4.3.1 to include non-identical transitions between ageing populations and a variable transition rate. This example shows how the state dependent distributed DDE framework addresses the inability of the linear chain technique to model dynamic ageing processes.

    Pérez-Ruixo et al. [34] studied the effect of recombinant human erythropoietin (EPO) on red blood cell precursors using a mathematical model. EPO is the protein responsible for controlling production of red blood cells and their precursors. The model arises from pharmacokinetic and pharmacodynamic data from patients receiving one dose of exogenous EPO. EPO was modelled through an open two compartment model of exogenous dose absorption and homoeostatic endogenous production rate, k_{EPO} , and the blood serum level (BSL). The bioavailable exogenous EPO was modelled as a dose dependent hyperbolic function satisfying

    \begin{equation*} F = F_0+ \frac{E_{max} {\rm Dose}}{ED_{50}+ {\rm Dose}}, \end{equation*}

    where {\rm Dose} is the amount of EPO administered. Exogenous EPO was absorbed through a dual absorption model into the depot and central compartments. The duration of first order absorption into the depot and central compartments are given by D_1 and D_2 , respectively. A fraction of the bioavailable exogenous EPO, f_r , was absorbed into the depot compartment before entering the central compartment at rate k_a . The depot concentration of EPO follows

    \begin{equation} \frac{ {\rm d}}{ {\rm dt}} A_1(t) = \left \{ \begin{array}{lll} \frac{ {\rm Dose} f_r F}{D_1} - k_a A_1 & {\rm if} & t \leqslant D_1\\ -k_a A_1 & {\rm if} & t \gt D_1, \end{array} \right. \end{equation} (5.1)

    The remaining exogenous EPO, (1-f_r)F enters the central compartment following a lag time t_{ {\rm lag}2} and is cleared linearly at the rate k_{20} . The volume of the central compartment is V_1 . The dynamics of exogenous EPO in the central compartment are given by

    \begin{equation} \frac{ {\rm d}}{ {\rm dt}} A_2(t) = \left \{ \begin{array}{lll} \frac{ {\rm Dose}(1- f_r) F}{D_2} + k_a A_1(t) + k_{32}A_3(t)- k_{23}A_2(t) & & \\ \quad - k_{20}A_2(t) + k_{epo} - \frac{V_{max}A_2(t)/V_2}{K_M+A_2(t)/V_2} & {\rm if} & t_{ {\rm lag}2} \leqslant t \leqslant D_2 \\ k_{epo} - \frac{V_{max}A_2(t)/V_1}{K_M+A_2(t)/V_1} & {\rm if} & t \gt D_2, t \lt t_{ {\rm lag}2}, \end{array} \right. \end{equation} (5.2)

    Finally, EPO enters the peripheral compartment from -and returns to- the central compartment linearly, so

    \begin{equation} \frac{ {\rm d}}{ {\rm dt}} A_3(t) = k_{23}A_2(t)-k_{32}A_3(t). \end{equation} (5.3)

    The total bioavailable EPO is given by

    \begin{equation*} C(t) = BSL + A_2(t)/V_1. \end{equation*}

    Pérez-Ruixo[34] considered 4 different pharmacodynamics models of erythrocyte response to exogenous EPO (titled the A, B, C and D models). In each of the 4 different pharmacodynamic models, the EPO dynamics are unchanged and described by equations (5.1), (5.2) and (5.3).

    Here, we describe the "B" model from [34]. Model B divides the erythrocyte progenitors, P(t) , into N_P compartments further subdivided into two distinct populations; EPO only affects the growth rate of the first population. Thus, the first N_P/2 compartments constitute the EPO sensitive population. Progression through these N_P compartments represents the ageing process of the progenitor cells. Once erythrocyte progenitors have reached maturity, they progress into the reticulocyte population. Once again, the maturation process of reticulocytes is modelled through a series of N_R transit compartments that are not sensitive to EPO. In this manner, the Pérez-Ruixo [34] model uses a concatenation of transit compartments to model the separate ageing processes of reticulocytes.

    The Pérez-Ruixo B model of erythrocyte progenitor and reticulocyte production is

    \begin{equation} \left. \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} P_1(t) & = k_{in}- \frac{S_{max}C(t)}{SC_{50}+C(t)}\frac{N_P}{T_P}P_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} P_i(t) & = \frac{S_{max}C(t)}{SC_{50}+C(t)} \frac{N_P}{T_P}\left[ P_{i-1}(t)-P_i(t)\right] \quad {\rm for} \quad i = 2, 3, ..., N_P/2 \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} P_{N_P/2+1}(t) & = \frac{S_{max}C(t)}{SC_{50}+C(t)}\frac{N_P}{T_P} P_{N_P/2}(t) - \frac{N_P}{T_P} P_{N_P/2+1}(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} P_i(t) & = \frac{N_P}{T_P}\left[ P_{i-1}(t)-P_i(t)\right] \quad {\rm for} \quad i = N_P/2+2, ..., N_P. \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} R_1(t) & = \frac{N_P}{T_P}P_{N_P}(t) - \frac{N_R}{T_R}R_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} R_i(t) & = \frac{N_R}{T_R}\left[ R_{i-1}(t)-R_i(t)\right] \quad {\rm for} \quad i = 2, 3, ..., N_R. \end{aligned} \right \} \end{equation} (5.4)

    By identifying the ingredients necessary from Remark 4.5, we will show how the distributed DDE framework from Section 4.3.1 can account for these separate ageing processes with distinct ageing velocities. Accounting for multiple ageing processes is not possible by rescaling time so the approach in [12] cannot be generalized to this case.

    The most immature erythrocyte progenitors are modelled by P_1(t) and are created from multipotent progenitors differentiating into the erythrocyte lineage at a constant rate k_{in} . Transit between the first N_P/2 compartments occurs at the variable rate

    \begin{equation*} V_e(t) = \frac{S_{max}C(t)}{SC_{50}+C(t)} \frac{N_P}{T_P} \quad {\rm with} \quad V_e^* = \frac{S_{max}BSL}{SC_{50}+BSL} \frac{N_P}{T_P}. \end{equation*}

    Using (4.22), we define \hat{V}_e(t) = V_e(t)/V_e^* , so the birth rate of precursor cells into P_2(t) is

    \begin{equation*} V_e(t) P_1(t) = \frac{\beta_e(t)}{\hat{V}_{e}(t)}. \end{equation*}

    Further, we see that the only removal of cells from the compartment model is due to transition to later compartments. Therefore, \mu(t) = 0 , and we have identified all the ingredients necessary in Remark 4.5. Therefore, for i = 2, 3, ...N_P/2 ,

    \begin{equation} P_i(t) = \int_{-\infty}^t \frac{V_e( \varphi)}{V_e^*} \ P_1( \varphi)g_{V_e^*}^{i}\left[\int_{ \varphi}^t \hat{V}_e(s) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation} (5.5)

    The N_P/2+1 st compartment satisfies

    \begin{equation*} \frac{ {\rm d}}{ {\rm dt}} P_{N_P/2+1}(t) = V_e(t) P_{N_P/2}(t) - \frac{N_P}{T_P} P_{N_P/2+1}(t). \end{equation*}

    Erythrocyte progenitors enter the first non-EPO sensitive ageing compartment, P_{N_P/2+1}(t) , with appearance rate

    \begin{equation*} \frac{\tilde{\beta}_e(t)}{V_p(t)} = V_e(t)P_{N/2}(t), \end{equation*}

    and then progress through the remaining N_P/2 compartments at a constant rate V_p(t) = V_p^* = N_P/T_P . Once again, we note that there is no removal of cells in any of the N_P/2 compartments, so \mu(t) = 0 . Further, since the ageing velocity is constant, \hat{V}_p^* = 1 . Therefore, a simple application of Remark 4.5 for constant ageing velocity, and using (5.5) gives

    \begin{align} P_{N_P}(t) & = \int_0^{\infty} \frac{\tilde{\beta}_e(t-\theta)}{N_P/T_P}g_{N_P/T_P}^{N_P/2}(\theta) \mathrm{d} \theta = \int_{-\infty}^{t} \frac{\tilde{\beta}_e(\theta)}{N_P/T_P}g_{N_P/T_P}^{N_P/2}(t-\theta) \mathrm{d} \theta \\[0.25cm] & = \int_{-\infty}^{t} \left[\frac{V_e(\theta)}{N_P/T_P} \int_{-\infty}^{\theta} V_e( \varphi) P_1( \varphi)g_{V_e^*}^{i-1}\left( \int_{ \varphi}^{\theta} \hat{V}_e(s) \mathrm{d} s\right) \mathrm{d} \varphi \right]g_{N_P/T_P}^{N_P/2}(t-\theta) \mathrm{d} \theta. \end{align} (5.6)

    Mature erythrocyte precursors enter into the most immature reticulocyte compartment, R_1(t) . Given (5.6), the differential equation for R_1(t) becomes

    \begin{equation*} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} R_1(t) & = \frac{N_P}{T_P}\underbrace{ \int_{-\infty}^{t} \left[\frac{V_e(\theta)}{N_P/T_P} \int_{-\infty}^{\theta} V_e( \varphi) P_1( \varphi)g_{V_e^*}^{i-1}\left( \int_{ \varphi}^{\theta} \hat{V}_e(s) \mathrm{d} s\right) \mathrm{d} \varphi \right]g_{N_P/T_P}^{N_P/2}(t-\theta) \mathrm{d} \theta}_{P_{N_P}(t)} \\ & \qquad {} - \frac{N_R}{T_R}R_1. \end{aligned} \end{equation*}

    Hence, the Pérez-Ruixo B model of reticulocyte production is equivalent to

    \begin{equation*} \begin{aligned} C(t) & = BSL + A_2(t)/V_1 \\ \frac{ {\rm d}}{ {\rm dt}} P_1(t) & = k_{in}- \frac{S_{max}C(t)}{SC_{50}+C(t)}\frac{N_P}{T_P}P_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} R_1(t) & = \frac{N_P}{T_P}\int_{-\infty}^{t} \left[\frac{V_e(\theta)}{N_P/T_P} \int_{-\infty}^{\theta} V_e( \varphi) P_1( \varphi)g_{V_e^*}^{N_p/2}\left( \int_{ \varphi}^{\theta} \hat{V}_e(s) \mathrm{d} s\right) \mathrm{d} \varphi \right]g_{N_P/T_P}^{N_P/2}(t-\theta) \mathrm{d} \theta \\ & \qquad {} - \frac{N_R}{T_R}R_1 \\ \frac{ {\rm d}}{ {\rm dt}} R_i(t) & = \frac{N_R}{T_R}\left[ R_{i-1}(t)-R_i(t)\right] \quad {\rm for} \quad i = 2, 3, ...N_R. \end{aligned} \end{equation*}

    Finally, we can use Remark 4.5 with the constant ageing velocity V_r(t) = V_r^* = N_R/T_R to solve the transit compartment system for R_i(t) to find

    \begin{equation} R_{i}(t) = \int_{0}^{\infty} \frac{T_R}{N_R} \beta_R(\sigma)g_{N_R/T_R}^{i}(\sigma) \mathrm{d} \sigma, \end{equation} (5.7)

    where

    \begin{equation*} \beta_R(\sigma) = \frac{N_P}{T_P}\int_{-\infty}^{\sigma} \left[\frac{V_e(\theta)}{N_P/T_P} \int_{-\infty}^{\theta} V_e( \varphi) P_1( \varphi)g_{V_e^*}^{N_p/2}\left( \int_{ \varphi}^{\theta} \hat{V}_e(s) \mathrm{d} s\right) \mathrm{d} \varphi \right]g_{N_P/T_P}^{N_P/2}(t-\theta) \mathrm{d} \theta. \\ \end{equation*}

    Using the techniques developed in Section 4.3.1, we have transformed the differential equations for the transit compartments for the erythrocyte progenitors and the reticulocytes into renewal type equations given by (5.6) and (5.7) [35]. Since Pérez-Ruixo [34] did not model reticulocyte mediated clearance of EPO, the cytokine and early progenitor dynamics are independent of the P_{N_P}(t) and R_{N_R}(t) concentrations. Consequently, the dynamics of equation (5.4) are completely determined by the dynamics of

    \begin{equation*} \begin{aligned} C(t) & = BSL + A_2(t)/V_1 \\ \frac{ {\rm d}}{ {\rm dt}} P_1(t) & = k_{in}- \frac{S_{max}C(t)}{SC_{50}+C(t)}\frac{N_P}{T_P}P_1(t), \end{aligned} \end{equation*}

    and the EPO concentrations given by equations (5.1), (5.2), and (5.3). We are now able to completely characterise the homeostatic behaviour of erythropoiesis by studying

    \begin{equation} \left. \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} A_1(t) & = -k_aA_1(t) \\ \frac{ {\rm d}}{ {\rm dt}} A_2(t) & = k_{epo} - \frac{V_{max}A_2/V_1}{K_M+A_2/V_1} \\ \frac{ {\rm d}}{ {\rm dt}} A_3(t) & = k_{23}A_2(t)-k_{32}A_3(t) \\ \frac{ {\rm d}}{ {\rm dt}} P_1(t) & = k_{in}- \frac{S_{max}C(t)}{SC_{50}+C(t)}\frac{N_P}{T_P}P_1(t), \end{aligned} \right \} \end{equation} (5.8)

    To ensure that the initial value problem (5.8) is equivalent to the Pérez-Ruixo model [34], we re-use the initial conditions for A_1(0), A_2(0), and A_3(0) . Since \mu = 0 and the initial conditions P_1(0) = P_i(0) are constant, we can set the history function for the progenitors, \rho_p(s) , to be \rho_p(s) = P_1(0) . The same can be done for the reticulocytes with \rho_r(s) = R_1(0) .

    We find the homeostatic concentration of EPO in the depot, central and peripheral compartments by solving

    \begin{equation*} \frac{ {\rm d}}{ {\rm dt}} A_1(t) = 0, \quad \frac{ {\rm d}}{ {\rm dt}} A_2(t) = 0, \quad \frac{ {\rm d}}{ {\rm dt}} A_3(t) = 0, \quad {\rm and} \quad \frac{ {\rm d}}{ {\rm dt}} P_1(t) = 0. \end{equation*}

    This yields the following homeostatic EPO concentrations

    \begin{equation*} A_1^* = 0, \quad A_2^* = \frac{V_1k_{epo}k_M}{V_{max}-k_{epo}}, \quad A_3^* = \frac{k_{23}}{k_{32}}A_2^*, \quad {\rm and} \quad C^* = BSL+A_2^*, \end{equation*}

    while the homeostatic progenitor concentration is

    \begin{equation*} P_1^* = \frac{k_{in}(SC_{50}+C^*)}{S_{max}C^*}\frac{T_P}{N_P}. \end{equation*}

    The simplified erythropoiesis dynamics (5.8) and homeostatic concentrations lead to the following proposition:

    Proposition 5.1. For positive parameter values, the homeostatic equilibrium point of equation (5.4) is locally asymptotically stable.

    Proof. The linearisation matrix of equation (5.8) about the equilibrium x^* = (A_1^*, A_2^*, A_3^*, P_1^*) is

    \begin{equation*} \mathbb{J}(x^*) = \left[ \begin{array}{cccc} -k_a & 0 & 0 & 0 \\ 0 & \frac{-V_{max}/V_1 k_M}{(k_M+A_2^*/V_1)^2} & 0 & 0 \\ 0 & k_{23} & -k_{32} & 0 \\ 0 & \frac{1}{V_1}\frac{S_{max}C^*}{(SC_{50}+C^*)^2} & 0 & -\frac{S_{max}C^*}{SC_{50}+C^*}\frac{N_P}{T_P}\\ \end{array} \right]. \end{equation*}

    The matrix \mathbb{J}(x^*) is lower triangular with strictly negative diagonal entries, so the eigenvalues are strictly negative and the equilibrium is locally asymptotically stable.

    This example illustrates how Remark 4.5 can be adapted to include a series of concatenated ageing processes. In the age structured PDE interpretation, each ageing process corresponds to a unique random variable modelling the transition between distinct stages. As we do not a priori expect the transition ages to be independent, interpreting the resulting ageing processes requires some care. The final renewal equation (5.8) includes a joint multivariate distribution representing the concatenation of distinct ageing processes.

    Further, Pérez-Ruixo [34] did not show that the homeostatic equilibrium is locally asymptotically stable. For the ODE system (5.4), the Jacobian would be a (3+N_P+N_R) \times (3+N_P+N_R) matrix with a degree (3+N_P+N_R) characteristic polynominal. In general, analytically finding the roots of a large degree polynominal is difficult. Hence, while the ODE (5.4) is obviously finite dimensional, it is analytically intractable.

    Conversely, the equivalent renewal equation (5.8) is simple to analyse and a similar argument to Proposition 3.1 shows that solutions of the renewal equation (5.8) evolving from non-negative initial conditions remain non-negative. The "A", "C" and "D" models can also be modelled as renewal equations through a simple application of the classical linear chain technique and the technique shown here.

    Roskos et al. [36] modelled the impact of exogenous administration of granulocyte colony stimulating factor (G-CSF) on neutrophil proliferation and maturation speed. G-CSF is a proinflammatory cytokine that binds to G-CSF specific receptors on mature neutrophil cells and controls neutrophil kinetics through a negative feedback loop [37,38]. G-CSF governs neutrophil production by increasing the effective proliferation of neutrophil precursors, reducing the maturation time of non-mitotic neutrophil precursors, and increasing release of neutrophil cells from the bone marrow into the blood. The dynamics of neutrophil production have been well-studied from both a mathematical and a pharmacometric point of view [5,12,24]. These models have used different techniques to incorporate the delays intrinsic to the system, such as discrete DDEs or transit compartment ODEs. Roskos et al. [36] model distinct stages of granulocyte production such as the bone marrow concentrations of metamyelocytes, M(t) ; band cells, B(t) ; and segmented neutrophil cells, S(t) . The ageing and maturation processes for each of these cell types is modelled through a series of three transit chains with N_M, N_B and N_S compartments, respectively. Moreover, band and segmented neutrophil cells can be shunted into circulation following the administration of G-CSF. We denote the metamyelocyte, band and segmented neutrophil cell shunting rates as \mu_m(t), \mu_b(t) and \mu_s(t)

    Administration of G-CSF is modelled in a similar way to the EPO model of Section 5.1 using a first order delayed absorption model. However, Roskos et al. [36] do not give the differential equations for exogenous administration of G-CSF other than to state that the clearance of G-CSF includes neutrophil receptor mediated clearance through the term

    \begin{equation*} CL_N/F = \frac{k_{cat}/F (B_p(t)+S_p(t))}{K_M+C(t)}, \end{equation*}

    where B_p(t) and S_p(t) are the number of circulating band and segmented neutrophil cells, respectively. Due to the feedback between the circulating neutrophil precursors and the cytokine C(t) , we are unable to completely reduce the Roskos model to a renewal type equation as was done in Section 5.1.

    The Roskos model for granulocyte production is

    \begin{align*} \frac{ {\rm d}}{ {\rm dt}} M_1(t) & = S_0+\frac{E_{mit}C(t)}{EC_{50}+C(t)} - \frac{N_M}{\tau_{meta}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)}\right) }M_1(t) \\ \frac{ {\rm d}}{ {\rm dt}} M_i(t) & = \frac{N_M}{\tau_{meta}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t) } \right) }\left( M_{i-1}(t)-M_{i}(t) \right) \quad {\rm for} \quad i = 2, ..., N_M \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} B_1(t) & = \frac{N_M}{\tau_{meta}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right)}M_{N_M}(t) - \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) } B_1(t) \\[0.2cm] & \qquad {} - \frac{E_{band}C(t)}{EC_{50}+C(t)}B_1(t) \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} B_i(t) & = \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) }\left[ B_{i-1}(t)-B_i(t)\right] - \frac{E_{band}C(t)}{EC_{50}+C(t)}B_i(t); \quad i = 2, ...N_B \\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} B_p(t) & = \sum\limits_{i = 1}^{N_B} \frac{E_{band}C(t)}{EC_{50}+C(t)}B_i(t)- (k_{\lambda}+ k_{bpmat})B_p(t) \\ \frac{ {\rm d}}{ {\rm dt}} S_1(t) & = \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right)}B_{N_B}(t) - \left( \frac{N_S}{\tau_{seg}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) } + \frac{E_{seg}C(t)}{EC_{50}+C(t)}\right) S_1(t)\\[0.2cm] \frac{ {\rm d}}{ {\rm dt}} S_i(t) & = \frac{N_S}{\tau_{seg}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) } \left[ S_{i-1}(t) - S_i(t) \right] - \frac{E_{seg}C(t)}{EC_{50}+C(t)}S_i(t); \quad i = 2, ..., N_S. \\ \frac{ {\rm d}}{ {\rm dt}} S_p(t) & = \sum\limits_{i = 1}^{N_S} \frac{E_{band}C(t)}{EC_{50}+C(t)}S_i(t)- (k_{\lambda}+ k_{bpmat})S_p(t), \end{align*}

    and is an example of a transit compartment model with variable ageing speed and linear clearance. The linear clearance terms are Hill type functions with a maximal clearance rate E_{j} given by

    \begin{equation*} \mu_j(t) = \frac{E_{j}C(t)}{EC_{50}+C(t)}. \end{equation*}

    Including these linear clearance terms in a transit compartment model is uncommon, but allows for the direct modelling of G-CSF mediated shunting of immature cells into circulation.

    By converting the model into a distributed DDE, we underline the link between clearance of cells in a transit compartment to the exponential decay present in the distributed DDE. Once again, we will proceed by identifying the ingredients discussed in Remark 4.5.

    As in Section 5.1, the most immature metamyelocytes ( M_1(t) ) are produced from the earlier progenitors at a constant baseline rate S_0 with the G-CSF dependent recruitment rate

    \begin{equation*} \frac{\beta_m(t)}{V_m(t)} = S_0 + \frac{E_{mit}C(t)}{EC_{50}+C(t)}. \end{equation*}

    Metamyelocytes progress through maturation at a G-CSF dependent rate

    \begin{equation*} V_m(t) = \frac{N_M}{\tau_{meta}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)}\right) }. \end{equation*}

    Metamyelocytes are not shunted into circulation following the administration of G-CSF, so \mu_m(t) = 0 . Therefore, the metamylocyte transit compartment model can be reduced to a distributed DDE using Remark 4.5 in an identical procedure to the Pérez-Ruixo model in Section 5.1. The most mature metamyelocyte population is given by

    \begin{equation} M_{N_M}(t) = \int_{-\infty}^t \frac{\beta_m(t)}{V_m(t)}g_{V_m^*}^{N_M}\left[\int_{ \varphi}^t \hat{V}_m(s) \mathrm{d} s\right] \mathrm{d} \varphi. \end{equation} (5.9)

    Immature neutrophil band cells, B_1(t) , are created at the birth rate

    \begin{equation*} \frac{\beta_b(t)}{\hat{V}_b(t)} = \frac{N_M}{\tau_{meta}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right)}M_{N_M}(t). \end{equation*}

    These band cells progress through the maturation compartments at the G-CSF dependent ageing rate

    \begin{equation*} V_b(t) = \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) } \quad {\rm with} \quad V_b^* = \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C^*}{EC_{50}+C^*} \right) }, \end{equation*}

    so the scaled ageing rate is \hat{V}_b(t) = V_b(t)/V_b^* . Inspecting the remaining terms in the equation for B_1(t) gives

    \begin{equation*} \mu_b(t) = \frac{E_{band}C(t)}{EC_{50}+C(t)}. \end{equation*}

    Therefore, using Remark 4.5, we find that the i -th band compartment satisfies

    \begin{equation} \begin{aligned} B_{i}(t) & = \int_{-\infty}^t \frac{\beta_b( \varphi)}{V_b( \varphi)} \exp\left[-\int_{ \varphi}^t \mu_b(s) \mathrm{d} s\right]g_{V_B^*}^{i}\left(\int_{ \varphi}^t \hat{V}_b(s) \mathrm{d} s\right) \mathrm{d} \varphi \end{aligned} \end{equation} (5.10)

    for i = 1, 2, ...N_B .

    Mature band cells, given by (5.10) with i = N_B , transition into the first segmented neutrophil cell compartment S_1(t) with creation rate

    \begin{equation*} \frac{\beta_s(t)}{\hat{V}_s(t)} = \frac{N_B}{\tau_{band}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right)}B_{N_B}(t) = V_b(t)B_{N_B}(t). \end{equation*}

    These cells transit through the segmented neutrophil population with G-CSF dependent ageing ( V_s(t) ) and clearance ( \mu_s(t) ) rates

    \begin{equation*} V_s(t) = \frac{N_S}{\tau_{seg}\left( 1-\frac{f_{mmt}C(t)}{EC_{50}+C(t)} \right) } \quad {\rm and} \quad \mu_s(t) = \frac{E_{seg}C(t)}{EC_{50}+C(t)}. \end{equation*}

    Therefore, we have identified all the ingredients in Remark 4.5 for the segmented neutrophil precursors, S(t) . The first segmented neutrophil cell compartment satisfies

    \begin{equation*} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} S_1(t) & = \overbrace{ V_b(t) \int_{-\infty}^t \frac{\beta_b( \varphi)}{V_b( \varphi)} \exp\left[-\int_{ \varphi}^t \mu_b(s) \mathrm{d} s\right]g_{V_B^*}^{N_b}\left(\int_{ \varphi}^t \hat{V}_b(s) \mathrm{d} s\right) \mathrm{d} \varphi}^{\beta_s(t)/\hat{V}_s(t)} \\ & \qquad {} -V_s(t) S_1(t) - \mu_s(t) S_1(t). \end{aligned} \end{equation*}

    Therefore, it is possible to replace the transit compartment system of ODEs for S_i(t) using Remark 4.5 to find

    \begin{align} S_i(t) & = \int_{-\infty}^t \overbrace{ V_b(\theta)\left[ \int_{-\infty}^{\theta} \frac{\beta_b( \varphi)}{V_b( \varphi)} \exp\left[-\int_{ \varphi}^{\theta} \mu_b(s) \mathrm{d} s\right] g_{V_B^*}^{N_b}\left(\int_{ \varphi}^{\theta} \hat{V}_b(s) \mathrm{d} s\right) \mathrm{d} \varphi\right]}^{\beta_s(\theta)/\hat{V}_s(\theta)} \\ & \qquad {} \times \exp\left[- \int_{\theta}^{t} \mu_s(x) \mathrm{d} x \right]g_{V_s^*}^{i}\left(\int_{\theta}^{t} \hat{V}_s(s) \mathrm{d} s\right) \mathrm{d} \theta \qquad {\rm for} \quad i = 1, 2, ..., N_s. \end{align} (5.11)

    The initial value problem studied in [36] was equipped with initial conditions for the cytokine equations as well as the N_M+N_S+N_B+2 compartments. Since \mu \neq 0 in general, to create an equivalent renewal type equation, we use the same initial conditions as [36] for the cytokine differential equations and follow [6] to construct appropriate history functions for M(t), B(t) and S(t) .

    Therefore, we can reduce the ODE model of granulopoiesis to a renewal-type equation with unchanged cytokine dynamics from [36] using the resulting DDEs for B_p(t) and S_p(t) . The resulting renewal equation is given by the equations describing the cytokine dynamics and the system of distributed DDEs

    \begin{equation*} \begin{aligned} \frac{ {\rm d}}{ {\rm dt}} B_p(t) & = \sum\limits_{i = 1}^{N_B} \frac{E_{band}C(t)}{EC_{50}+C(t)}B_i(t)- (k_{\lambda}+ k_{bpmat})B_p(t) \\ \frac{ {\rm d}}{ {\rm dt}} S_p(t) & = \sum\limits_{i = 1}^{N_S} \frac{E_{band}C(t)}{EC_{50}+C(t)}S_i(t)- (k_{\lambda}+ k_{bpmat})S_p(t), \end{aligned} \end{equation*}

    where B_i(t) and S_i(t) are given by (5.10) and (5.11), respectively.

    In this example, we have shown how to concatenate multiple ageing processes with distinct ageing velocities, as well as how to include the loss of cells throughout the ageing process. Once again, we can use a similar argument to Proposition 3.1 to ensure that the solutions evolving from non-negative initial data remain non-negative.

    In this work, we have shown how to reduce age structured PDEs to possibly state-dependent DDEs. Our derivation shows how the correction factor discussed in Section 2.1 results naturally from considering the hazard rate at which cells exit maturation, and generalises the derivation of Craig et al. [5] to the non-deterministic case.

    In Section 3, we analysed the general distributed DDE that arises from the age structured population model. We showed, in Proposition 3.1, that populations evolving from non-negative initial conditions remain non-negative, regardless of the density K_A(t) . By linearising the distributed DDE, we showed, in Proposition 3.2, that stability analysis of the general DDE is analytically tractable. We characterized the stability of a generic equilibrium solution as a function of the linearisation of the growth function F(x^*, \bar{x}^*) .

    Next, we considered the state-dependent DDE in the case of the degenerate, uniform and gamma distributions. Choosing a degenerate distribution leads to the familiar state-dependent discrete DDE, while uniformly distributed DDEs are reducible to discrete DDEs with two state dependent delays. Finally, in the case of gamma distributed DDEs, we explicitly related transit compartment models that include variable transit rates with gamma distributed DDEs in Theorem 4.4. As shown in [12], it can be simpler to analyse stability of equilibria and positivity of solutions of a distributed DDE than the corresponding ODE. However, the ODE models may be simpler to simulate numerically. The equivalence between the differential equations allows for the resulting model to be analysed in the more convenient setting.

    By the means of two examples, we showed how to express transit compartment models as an equivalent DDE or renewal equation. First, we showed how to incorporate a variable transit rate into a distributed DDE using a simple application of Theorem 4.4. Next, we demonstrated that our method is capable of including multiple distinct ageing processes in the form of a multivariate distributed DDE. Lastly, we showed how a linear clearance term in each of the transit compartments can be included in the equivalent DDE model. Analysis of the renewal equation was shown to be simpler than the corresponding ODE system, and we were able to easily characterise the stability of the homeostatic equilibria.

    This work emphasizes the link between transit compartment ODEs and delay differential equations. While this link has been known for over 50 years, we explicitly establish it for compartment models with variable transit rates. We demonstrated that these transit compartment models are equivalent to state dependent distributed DDEs. The equivalence between easy-to-simulate ODE models and the simpler to analyse distributed DDEs allows modellers to use the formulation that is most convenient for their purposes. Consequently, the framework developed in this article allows for researchers to incorporate both external control of ageing rates and heterogeneous, non-deterministic maturation age into models of physiological maturation processes.

    TC would like to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) for funding through the PGS-D program and the Alberta Government for funding through the Sir James Lougheed award of distinction. MC and ARH are grateful for funding through the NSERC Discovery Grant program.

    All authors declare no conflicts of interest in this paper.



    [1] A. G. McKendrick, Applications of mathematics to medical problems, Proc. Edinburgh Math. Soc., 44 (1925), 98–130.
    [2] E. Trucco, Mathematical models for cellular systems. The Von Foerster equation. Part II, Bull. Math. Biophys., 27 (1965), 449–471.
    [3] J. A. J Metz and O. Diekmann, The Dynamics of Physiologically Structured Populations, vol. 68 of Lecture Notes in Biomathematics. Berlin, Heidelberg: Springer, 3 ed., 1986.
    [4] H. L. Smith, Reduction of structured population models to threshold-type delay equations and functional differential equations: a case study, Math. Biosci., 113 (1993), 1–23.
    [5] M. Craig, A. R. Humphries and M. C. Mackey, A mathematical model of granulopoiesis incor-porating the negative feedback dynamics and kinetics of G-CSF/neutrophil binding and internal-ization, Bull. Math. Biol., 78 (2016), 2304–2357.
    [6] T. Cassidy and A. R. Humphries, A mathematical model of viral oncology as an immuno-oncology instigator, Math. Med. Biol., Online First (2019), dqz008. Available from: https://doi.org/10.1093/imammb/dqz008.
    [7] A. Otto and G. Radons, Transformations from variable delays to constant delays with applica-tions in engineering and biology, in Time Delay Syst. (T. Insperger, T. Ersal, and G. Orosz, eds.), vol. 7 of Advances in Delays and Dynamics, 169–183, Cham: Springer International Publishing, 2017.
    [8] S. Bernard, Moving the boundaries of granulopoiesis modelling, Bull. Math. Biol., 78 (2016), 2358–2363.
    [9] J. M. Mahaffy, J. Bélair and M. C. Mackey, Hematopoietic model with moving boundary condi-tion and state dependent delay: application in erythropoiesis, J. Theor. Biol., 190 (1998), 135–146.
    [10] D. R. Cox, Regression models and life-tables, J. R. Stat. Soc., 34 (1972), 187–220.
    [11] E. L. Kaplan and P. Meier, Nonparametric estimation from incomplete observations, J. Am. Stat. Assoc., 53 (1958), 457.
    [12] D. Câmara de Souza, M. Craig, T. Cassidy, et al., Transit and lifespan in neutrophil production: implications for drug intervention, J. Pharmacokinet. Pharmacodyn., 45 (2018),59–77.
    [13] J. K. Hale and S. M. Verduyn Lunel, Introduction to Functional Differential Equations, vol. 99 of Applied Mathematical Sciences. New York, NY: Springer New York, 1993.
    [14] Y. Hino, S. Murakami and T. Naito, Functional Differential Equations with Infinite Delay, vol. 1473 of Lecture Notes in Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg, 1991.
    [15] W. Liu, T. Hillen and H. I. Freedman, A mathematical model for M-phase specific chemotherapy including the G0-phase and immunoresponse., Math. Biosci. Eng., 4 (2007), 239–259.
    [16] F. Hartung, T. Krisztin, H. O. Walther, et al., Chapter 5 Functional Differential Equations with State-Dependent Delays: Theory and Applications, in Handb. Differ. Equations (A. Canada, P. Drabek, and A. Fonda, eds.), ch. 5, pp. 435–545, North Holland 2004: Elsevier, 3rd ed., 2006.
    [17] Y. Yuan and J. Bélair, Stability and Hopf Bifurcation Analysis for Functional Differential Equa-tion with Distributed Delay, SIAM J. Appl. Dyn. Syst., 10 (2011), 551–581.
    [18] H. L. Smith, An Introduction to Delay Differential Equations with Applications to the Life Sci-ences, vol. 57 of Texts in Applied Mathematics. New York, NY: Springer New York, 2011.
    [19] A. Teslya, Predator-Prey Models With Distributed Time Delay. Doctoral dissertation, McMaster University, 2015.
    [20] T. Vogel, Systèmes Déferlants, Systèmes Héréditaires, Systèmes Dynamiques, in Proc. Int. Symp. Nonlinear Vib., (Kiev), pp. 123–130, Academy of Sciences USSR, 1961.
    [21] N. MacDonald, Time Lags in Biological Models. Berlin: Springer, 1978.
    [22] W. Krzyzanski, Interpretation of transit compartments pharmacodynamic models as lifespan based indirect response models., J. Pharmacokinet. Pharmacodyn., 38 (2011),179–204.
    [23] W. Gurney, R. Nisbet and S. Blythe, The systematic formulation of models of stage-structured populations, in Dyn. Physiol. Struct. Popul. (J. A. J. Metz and O. Diekmann, eds.), ch. 11, 474–493, Berlin, Heidelberg: Springer Berlin Heidelberg, 3 ed., 1986.
    [24] A. L. Quartino, M. O. Karlsson, H. Lindman, et al., Characterization of endogenous G-CSF and the inverse correlation to chemotherapy-induced neutropenia in patients with breast cancer using population modeling, Pharm. Res., 31 (2014), 3390–3403.
    [25] L. Glass, Dynamical disease: Challenges for nonlinear dynamics and medicine, Chaos An Inter-discip. J. Nonlinear Sci., 25 (2015).
    [26] S. Rubinow and J. Lebowitz, A mathematical model of neutrophil production and control in normal man, J. Math. Biol., 225 (1975), 187–225.
    [27] M. C. Mackey, Unified hypothesis for the origin of aplastic anemia and periodic hematopoiesis, Blood, 51 (1978), 941–956.
    [28] C. Colijn and M. C. Mackey, A mathematical model of hematopoiesis - I. Periodic chronic myelogenous leukemia, J. Theor. Biol., 237 (2005), 117–132.
    [29] F. Crauste and M. Adimy, Modeling and asymptotic stability of a growth factor-dependent stem cell dynamics model with distributed delay, Discret. Contin. Dyn. Syst. - Ser. B, 8 (2007), 19–38.
    [30] T. Hearn, C. Haurie and M. C. Mackey, Cyclical neutropenia and the peripherial control of white blood cell production, J. Theor. Biol., 192 (1998), 167–181.
    [31] L. E. Friberg, A. Henningsson, H. Maas, et al., Model of chemotherapy-induced myelosuppres-sion with parameter consistency across drugs, J. Clin. Oncol., 20 (2002), 4713–4721.
    [32] G. von Schulthess and N. Mazer, Cyclic neutropenia (CN): A clue to the control of granu-lopoiesis, Blood, 59 (1982), 27–37.
    [33] W. Krzyzanski, P. Wiczling, P. Lowe, et al., Population modeling of filgrastim PK-PD in healthy adults following intravenous and subcutaneous administrations, J. Clin. Pharmacol., 50 (2010), 101S–112S.
    [34] J. J. Pérez-Ruixo, W. Krzyzanski and J. Hing, Pharmacodynamic analysis of recombinant human erythropoietin effect on reticulocyte production rate and age distribution in healthy subjects., Clin. Pharmacokinet., 47 (2008), 399–415.
    [35] O. Diekmann, M. Gyllenberg and J. A. J. Metz, Finite dimensional state representation of linear and nonlinear delay systems, J. Dyn. Differ. Equations, 30 (2018), 1439–1467.
    [36] L. K. Roskos, P. Lum, P. Lockbaum, et al., Pharmacokinetic/pharmacodynamic modeling of pegfilgrastim in healthy subjects, J. Clin. Pharmacol., 46 (2006), 747–757.
    [37] A. Roberts, G-CSF: A key regulator of neutrophil production, but that's not all!, Growth Factors, 23 (2005), 33–41.
    [38] E. Shochat, V. Rom-Kedar and L. Segel, G-CSF control of neutrophils dynamics in the blood., Bull. Math. Biol., 69 (2007), 299–338.
  • This article has been cited by:

    1. Sean T. Vittadello, Scott W. McCue, Gency Gunasingh, Nikolas K. Haass, Matthew J. Simpson, A novel mathematical model of heterogeneous cell proliferation, 2021, 82, 0303-6812, 10.1007/s00285-021-01580-8
    2. Tyler Cassidy, Antony R. Humphries, Morgan Craig, Michael C. Mackey, Characterizing Chemotherapy-Induced Neutropenia and Monocytopenia Through Mathematical Modelling, 2020, 82, 0092-8240, 10.1007/s11538-020-00777-0
    3. Tyler Cassidy, Morgan Craig, Aaron Goldman, Determinants of combination GM-CSF immunotherapy and oncolytic virotherapy success identified through in silico treatment personalization, 2019, 15, 1553-7358, e1007495, 10.1371/journal.pcbi.1007495
    4. Tyler Cassidy, Daniel Nichol, Mark Robertson-Tessi, Morgan Craig, Alexander R. A. Anderson, Dominik Wodarz, The role of memory in non-genetic inheritance and its impact on cancer treatment resistance, 2021, 17, 1553-7358, e1009348, 10.1371/journal.pcbi.1009348
    5. Jong Hyuk Byun, Yunil Roh, In-Soo Yoon, Kwang Su Kim, Il Hyo Jung, Antony R Humphries, Fractional transit compartment model for describing drug delayed response to tumors using Mittag-Leffler distribution on age-structured PKPD model, 2022, 17, 1932-6203, e0276654, 10.1371/journal.pone.0276654
    6. Tomáš Gedeon, Antony R. Humphries, Michael C. Mackey, Hans-Otto Walther, Zhao Wang, Operon dynamics with state dependent transcription and/or translation delays, 2022, 84, 0303-6812, 10.1007/s00285-021-01693-0
    7. Tyler Cassidy, Peter Gillich, Antony R Humphries, Christiaan H van Dorp, Numerical methods and hypoexponential approximations for gamma distributed delay differential equations, 2022, 87, 0272-4960, 1043, 10.1093/imamat/hxac027
    8. Tyler Cassidy, Distributed Delay Differential Equation Representations of Cyclic Differential Equations, 2021, 81, 0036-1399, 1742, 10.1137/20M1351606
    9. Zhao Wang, Antony R. Humphries, Numerical Methods for Delay Differential Equations with Threshold State-Dependent Delay, 2022, 55, 24058963, 163, 10.1016/j.ifacol.2022.11.351
    10. Yueli Huang, Jin-E Zhang, Asymptotic stability of impulsive stochastic switched system with double state-dependent delays and application to neural networks and neural network-based lecture skills assessment of normal students, 2024, 9, 2473-6988, 178, 10.3934/math.2024011
    11. Tyler Cassidy, A Continuation Technique for Maximum Likelihood Estimators in Biological Models, 2023, 85, 0092-8240, 10.1007/s11538-023-01200-0
    12. Justin Le Sauteur-Robitaille, Powel Crosley, Mary Hitt, Adrianne L. Jenner, Morgan Craig, Mathematical modeling predicts pathways to successful implementation of combination TRAIL-producing oncolytic virus and PAC-1 to treat granulosa cell tumors of the ovary, 2023, 24, 1538-4047, 10.1080/15384047.2023.2283926
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5843) PDF downloads(1036) Cited by(12)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog