Histological analysis of Xylella fastidiosa infection in Quercus pyrenaica in Northern Portugal

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

1.   Introduction

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 $K_A(t)$ through

which satisfies

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 $K_A(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 $K_A(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 $K_A(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.

2.   From McKendrick type equations to state dependent delays

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 $V_a(t)$, where $V_a(t)$ satisfies

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

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

The death rate of immature individuals is given by $\mu(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 $K_A(t)$. The transition rate, $h(a(t))$, is the instantaneous change in probability that an individual matures at age $a(t+\Delta t)$, given that the individual has not matured at age $a(t)$. Formally, using the definition of conditional probability,

Multiplying by unity gives

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

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

where the integral term

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

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)

with initial data

Proof. The characteristics of equation (2.1) satisfy

and hence are given by

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

Equation (2.7) has solution

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

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

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) = \rho(s)$ for $s < t_0$ for the DDE. To do this, it is natural to assume that an an immature individual with positive age $a > 0$ at time $t_0$ was born at sometime $s < t_0$. Since the PDE (2.1) is not defined for $s < t_0$, we are free to prescribe fixed values for $V_a(s) = V_a^*$ and $\mu(x(s)) = \mu^*$ for $s < t_0$. Then, imposing that individuals born at time $s < t_0$ evolved according to the McKendrick Equation, we have $a = V_a^*(t_0-s)$, or $s = t_0-a/V_a^*$. Hence, the initial condition $f(a)$ defines the history function $\rho$ through

Therefore defining $x(s) = \rho(s)$ this way, for $s < t_0$, 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 $T_0 = t- \varphi$ have age

for $\varphi > 0.$ So we have

At time $t$, the rate at which individuals mature is

By defining, for any density $K_A(t)$,

we have

Consequently, using (2.11) and defining the history $\rho(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).

2.1.   Accounting for the random maturation threshold

Further inspection of equation (2.10) reveals a ratio of ageing speeds $V_a(t)/V_a(t- \varphi)$ in the integral term

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 $V_a(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].

3.   Properties of state dependent delay differential equations

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

and consider the IVP

where $F(x, y) \in\mathcal{C}^1(\mathbb{R}^2, \mathbb{R})$ and $\gamma(x(t)) \in \mathcal{C}^1(\mathbb{R}, \mathbb{R})$ with

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

and $\lambda$ is the Lebesgue measure on $\mathbb{R}$. $L_1(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 $L_1(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 $L_1(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 $\gamma(x(t))$ satisfy equations (3.3). Moreover, assume that the history function satisfies

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

Proof. As $\rho(s) \geqslant 0$, it is simple to see that

We have a series of cases.

1) If $\rho(t_0) = x(t_0) > 0$ then $F(x(t_0), \bar{x}(t_0)) > 0$. Therefore,

and using Gronwall's inequality, we have

2) If $\rho(t_0) = 0$ and $\rho(s) = 0$ $A$-almost everywhere in $(-\infty, t_0)$, then $x(t) = 0$ is the solution of the IVP.

3) Finally, if $\rho(t_0) = 0$ and $\rho(s) > 0$ on a set of $A$-positive measure in $(-\infty, t_0)$ then $\bar{x}(t_0) > 0$ and

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 > t_0$.

3.1.   Linearisation of the DDE

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

At the equilibrium $x(t) = x^*$, $V_a(t) = V_a^*$, so

then, by evaluating (3.1) with $x(t) = x^*$, the homeostatic delayed term $\bar{x}^*$ in (3.4) satisfies

where $\mathcal{L}[f](s)$ is the Laplace transform of $f(x)$ evaluated at $s$.

Hence, $\bar{x}^*$ is a function of the density $K_A(t)$. However, if desired, it is possible to vary the homeostatic death rate $\mu^*$ to ensure that the equilibria value $x^*$ does not change for different densities $K_A(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 $V_a(t) = V_a^*$ and $\mu(s) = \mu^*$. By expanding exponential integral to leading order following ^[6], it is possible to relax the assumption that $\mu(s) = \mu^*$. Now, define $\bar{z}(t) = \bar{x}(t)- \bar{x}^*$ so that

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

By making the ansatz

we compute the expression for $\bar{z}(t)$ from (3.5)

Therefore

where $\gamma^* = \partial_x\gamma(x(t))|_{x = x^*}$, $k_1 = \partial_a F(a, b)|_{(x, \bar{x})}$ and $k_2 = \partial_{b} F(a, b)|_{(x, \bar{x})}$. Dropping nonlinear terms, the linearised equation is

The characteristic equation corresponding to (3.6) is

Through a standard analysis, we study the local stability of the equilibrium $x^*$ for a density $K_A(t)$.

Proposition 3.2. 1) If

the equilibrium point $x^*$ is locally asymptotically stable.

2) If

the equilibrium point $x^*$ is unstable.

Proof. 1) Let $\lambda^*$ be a root of (3.7) and assume for contradiction that $\Re(\lambda^*) \geqslant 0.$ We necessarily have

and we calculate

We note that

While, for arbitrary $\nu = \nu_r+i\nu_i \in \mathbb{C}$,

Moreover, if $\nu_r \geqslant 0$,

Therefore, using the assumption in 1), we find

which is a contradiction, so no such $\lambda^*$ 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

and note that $g(\lambda)$ is continuous with

Then, there must be a real $\lambda^* > 0$ such that $g(\lambda^*) = 0$. The equilibrium is therefore unstable.

We note that if $k_2 > 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 $k_2 < 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 $\mu = \gamma = 0$.

4.   Distributed delay differential equations with specific maturation probabilities

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 $K_A(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 $x^*\in L_1(A)$ is an equilibrium point so that $\mu(x^*) = \mu^*$ and $V_a(t) = V_a^*$. Denote the homeostatic maturation time as the first moment of the random variable $A$ with constant ageing rate $V_a^*$,

Consequently, the expected homeostatic maturation age is given by $\mathcal{T} = V_a^* \tau^*$.

We first consider the degenerate distribution concentrated at $\mathcal{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.

4.1.   Deterministic maturation

Assuming that maturation is a deterministic process and occurs after achieving the threshold age $\mathcal{T}$ implies that $K_A(t)$ is the degenerate distribution concentrated at $\mathcal{T}$ with

where $\delta(x)$ is the Dirac delta function. In the deterministic case, all individuals mature at precisely the same age $\mathcal{T}$. At the equilibrium $x^*$, using (3.7), the characteristic equation is

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 $K_A(t)$ given by (4.1), the threshold maturation age $\mathcal{T}$ allows us to calculate when an individual that matures at time $t$ began maturation. The maturation time, $\tau(x(t))$, must satisfy the implicit threshold condition

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

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

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

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

In the case that $\rho(s) = x^*$, then it is simple to set $\tau(0) = \tau^*$. However, for more general initial data $\rho(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].

4.2.   Uniformly Distributed Maturation

We consider uniformly distributed DDEs centered about the expected homeostatic maturation age $\mathcal{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 $\mathcal{T}$ is

At the equilibrium $x^*$, with $K_A(t)$ given by the uniform density (4.6), the characteristic equation (3.7) is

$\mathcal{T}-V_a^*\delta$ and $\mathcal{T}+V_a^*\delta$ 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, $\tau_{min}(x(t))$ and $\tau_{max}(x(t))$, are state dependent, and implicitly defined by

We note that, at homeostasis, $V_a(s) = V_a^*$ so

Recalling that $\mathcal{T} = V_a^*\tau^*$, the terms $\tau^* - \delta$ and $\tau^*+\delta$ 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

Thus the state dependent uniform distributed DDE is

4.2.1.   Reduction to Discrete DDE

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. $A_U(t)$ satisfies the differential equation

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

We note that

so that, integrating by parts,

Consequently, the derivative of $A_U(t)$ is

To finish the proof, we note that, similar to (4.5), $\tau_{min}(x(t))$ and $\tau_{max}(x(t))$ solve the following differential equations

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

By writing the delay term $A_U(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

with suitably chosen initial data.

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

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

and

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

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)$.

4.3.   Gamma distributed maturation and a generalized linear chain technique

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

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

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

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

Using a common denominator gives

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

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

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

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

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

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

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.

4.3.1.   A Generalized linear chain technique

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}$,

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

Proof. We first note that

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

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

Now we define,

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

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

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

Importantly, Lemma 4.3 ensures that

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

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,$

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

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

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,

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.

5.   Examples from hematopoiesis

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.

5.1.   Pérez-Ruixo model of reticulocyte production

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

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

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

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

The total bioavailable EPO is given by

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

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

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

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$,

The $N_P/2+1$st compartment satisfies

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

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

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

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

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

where

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

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

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

This yields the following homeostatic EPO concentrations

while the homeostatic progenitor concentration is

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

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.

5.2.   Roskos's model of granulocyte production

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

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

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

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

Metamyelocytes progress through maturation at a G-CSF dependent rate

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

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

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

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

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

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

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

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

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

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

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.

6.   Discussion and conclusion

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.

Acknowledgments

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.

Conflict of interest

All authors declare no conflicts of interest in this paper.