Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase
Volume 1, Number 4, December 2020, pp.412-424 https://doi.org/10.5206/mase/10857

A PROCEDURE FOR DERIVING NEW ODE MODELS: USING THE

GENERALIZED LINEAR CHAIN TRICK TO INCORPORATE PHASE-TYPE

DISTRIBUTED DELAY AND DWELL TIME ASSUMPTIONS.

PAUL J. HURTADO AND CAMERON RICHARDS

Abstract. Ordinary differential equations (ODE) models are used in a wide variety of applications

throughout the sciences. Despite their widespread use, these models are sometimes viewed as inflexible

with respect to time delay and dwell time assumptions. The Generalized Linear Chain Trick (GLCT)

provides a theoretical foundation that can help modelers incorporate much more flexible phase-type

distributed delay (or dwell time) assumptions into ODE models, including traditional exponential

and Erlang distribution assumptions. The GLCT serves as a bridge between stochastic processes

and dynamical systems theory for ODEs, opening up opportunities to use concepts and tools from

Markov chain theory in the development, analysis, and interpretation of ODE models. To facilitate

the practical application of this theory, in this paper we introduce a new GLCT-based procedure for

deriving new ODE models by generalizing or approximating existing ODE, DDE, or distributed delay

equation models. We apply this procedure to multiple models from the literature, using it to derive

new models that are generalizations or approximations of those models.

1. Introduction

Ordinary differential equations (ODE) models are widely used in the sciences, in part because of the

relative ease of formulating and analyzing ODE models [35, 25, 3, 7, 1]. However, they are often criticized

for their limited capacity to only incorporate a narrow range of delay and dwell time assumptions. Such

delays are more easily incorporated into models using other mathematical frameworks, e.g., delays

can be modeled most generally using integral equations or integro-differential equations to incorporate

distributed delays into dynamic models, or using delay differential equations (DDEs) to incorporate

fixed delays [37, 9, 8, 31, 19, 30]. In an ODE framework, the linear chain trick has long been used

to incorporate exponential and Erlang distributed delays [23, 22, 26, 34, 6]. However, recently this

technique has been generalized to a much broader family of delay or dwell time distributions [16].

This generalized linear chain trick (GLCT) theory [16] allows modelers to incorporate delays and

dwell times that obey phase-type distributions. The phase-type family of distributions is made up

of the set of all possible absorption time distributions for continuous time Markov chains (CTMCs)

with one or more transient states and at least one absorbing state. These include exponential, Erlang

(gamma with an integer shape parameter), hyperexponential, hypoexponential (generalized Erlang), and

Coxian distributions [4, 12, 27, 14]. The GLCT also permits the use of similar time-varying versions of

such distributions [16]. Together, these tools and techniques enable modelers to draw from a richer set

of ODE model assumptions when constructing new models, and a framework for more clearly seeing

Received by the editors 1 July 2020; accepted 15 December 2020; published online 29 December 2020.

2010 Mathematics Subject Classification. Primary 92B99; Secondary 37N25, 37M05.
Key words and phrases. linear chain trick; gamma chain trick; phase-type distribution; Coxian distribution; Hypoexpo-

nential distribution.
PJH was supported in part by the Sloan Scholars Mentoring Network of the Social Science Research Council with

funds provided by the Alfred P. Sloan Foundation; and NSF Grant No. DEB-1929522.

412

https://ojs.lib.uwo.ca/mase
https://dx.doi.org/https://doi.org/10.5206/mase/10857


A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 413

how underlying stochastic model assumptions are reflected in the structure of mean field ODE models.

Moreover, concepts from Markov chain theory, including reward process theory, can be used in the

analysis of ODEs derived using the GLCT as well as the interpretation of general analytical results [18].

To help modelers put the GLCT into practice, in this paper we introduce a new GLCT-based procedure

that can be used to quickly formulate new ODE models that explicitly incorporate phase-type distributed

delays, or that otherwise make use of Markov chains to organize complex state transition structures

in such models. Below, we describe this relatively simple procedure and illustrate it’s application by

using it to derive new ODE models that are extensions of models taken from the literature. These

include ODE models with no explicit delays, DDE models, and distributed delay models in the form of

integro-differential equations.

In section 1.1 we review the GLCT framework for phase-type distributions, and the standard Linear

Chain Trick (LCT). In section 1.2, we review procedures for deriving ODE models as special cases of

distributed delay equations, or as approximations of DDEs. In section 2.1 we describe the GLCT-based

procedure. We then use it to extend multiple models starting in section 2.2, where we generalize a

tumor growth inhibition model by Simeoni et al [32]. In section 2.3 we generalize an opioid epidemic

model by Battista et al [2], then a within-host immune-pathogen interaction model by Hurtado [15] in

section 2.4, and a model of cell-to-cell HIV spread by Culshaw et al [5] in section 2.5.

1.1. Generalized Linear Chain Trick. We here provide a statement of the Generalized Linear Chain

Trick (GLCT) for phase-type distributions. The GLCT in its most general form, as detailed in [16],

extends Lemma 1.1 below to include time-varying parameters (analogous to extending homogeneous

Poisson process rates to the time-varying rates of inhomogeneous Poisson processes).

Phase-type distributions are a family of continuous, matrix exponential distributions that can be

thought of as the absorption time distributions for CTMCs with one or more absorbing states. They

are parameterized by a column vector v, which is the initial distribution vector across the set of

transient states1, and the transient state block M of the transition rate matrix, which together define

the corresponding CTMC. Equations for the density function f(t), cumulative distribution function

F(t), and jth moments E(Y j) of a phase-type distributed random variable Y are given in eqs. (1.1)

below, where superscript T indicates transposition and 1 is an appropriately long column vector of ones

(for more on phase-type distributions, see [17, 18, 4, 14, 12, 27]).

f(t) = vT eMt (−M1) (1.1a)

F(t) = 1−vT eMt 1 (1.1b)

E(Y j) = j! vT (−M)−j1. (1.1c)

Lemma 1.1 (GLCT for phase-type distributions [16]). Assume individuals enter a state (call

it state X) at rate I(t) ∈ R and that the distribution of time spent in state X follows a continuous
phase-type distribution given by the length k initial probability vector v and the k ×k matrix M. Then
partitioning X into k sub-states Xi, and denoting the corresponding amount of individuals in state Xi at

time t by xi(t), then the mean field equations for the quantities xi are given by

d

dt
x(t) = vI(t) + MT x(t) (1.2)

where the rate of individuals leaving each of these sub-states of X is given by the vector (−M 1) ◦ x,
where ◦ is the Hadamard (element-wise) product of the two vectors, and thus the total rate of individuals
leaving state X is given by the sum of those terms, i.e., by (−M 1)Tx = −1TMTx.

1The full initial distribution vector, for a CTMC with one absorbing state, would be [vT , v0]
T .



414 P. HURTADO AND C. RICHARDS

The standard Linear Chain Trick (LCT) is well known (see [16] and references therein) and is a

special case of Lemma 1.1 above. However, it is usually stated without the above matrix-vector notation.

The following is a formal statement of the standard LCT, but here we have slightly generalized it to

include generalized Erlang distributions (i.e., the sum of k independent exponentially distributed random

variables, each with potentially different rates ri) as this only changes a few subscripts in the mean field

equations. See [16, 34] for a similar statement of the standard LCT (for Erlang distributions).

Corollary 1.1 (Linear Chain Trick for Hypoexponential Distributions). Consider the GLCT

above. Assume that the dwell-time distribution is a generalized Erlang (hypoexponential) distribution

with rates r = [r1, r2, . . . ,rk]
T, where ri > 0, or an Erlang distribution with rate r (all ri = r) and

shape k (or if written in terms of shape k and mean τ = k/r, use r = k/τ). Then the corresponding

mean field equations are

dx1
dt

= I(t)−r1 x1
dx2
dt

= r1 x1 −r2 x2
...

dxk
dt

= rk−1 xk−1 −rk xk.

(1.3)

Proof. The phase-type distribution formulation of this generalized Erlang distribution is given by

v =



1

0
...

0


 and M =




−r1 r1 0 · · · 0

0 −r2 r2
... 0

...
...

...
...

...

0 0
... −rk−1 rk−1

0 0 · · · 0 −rk



. (1.4)

Substituting these into eq. (1.2) yields eqs. (1.3). If all ri = r then the phase-type distribution is an

Erlang distribution with rate r and shape k (mean k/r and coefficient of variation 1/
√
k). �

1.2. Deriving ODEs from Delay Differential Equations & Distributed Delay Equations. The

first step in our GLCT-based procedure detailed below is to take an existing ODE, DDE, or distributed

delay equation and derive ODEs that either approximate (in the DDE case) or are a special case of (in

the distributed delay equation case) that model. Here we briefly review how this can be done.

First, for distributed delay equations (e.g., integral equations) the standard linear chain trick can be

applied, as detailed in [16]. This involves assuming Erlang distributed delays, and from those distributed

delay equations deriving an ODE model by differentiating integral equations and employing the recursive

property of Erlang density functions. Alternatively, in some cases you may be able to use results such

as Lemma 1.1 or Corollary 1.1 to formulate these equations directly.

For DDEs, one can use the following (well known) technique to obtain an approximating system

of ODEs. Here we summarize this approximation as described in section 7.3 of [34] (for more on this

technique, see [10, 11, 24]). Consider an equation of the form

d

dt
x(t) = f(x(t),x(t− τ)) (1.5)

with solution x(t) that includes the initial data over [−τ,0] given by x(t) = φ(t), t ∈ [−τ,0].
Our goal is to obtain a system of ODEs that approximates eq. (1.5). To do this, first observe that

the assumed discrete delay of τ time units is equivalent to assuming a distributed delay model where



A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 415

the delay distribution is a Dirac-δ function shifted by τ time units, δ(t− τ). That is, the CDF for this
distribution is F(t) = 0 for t < τ and then it jumps to F(t) = 1 for t ≥ τ.

This Dirac-δ distributional assumption can be approximated with an Erlang (gamma) distribution

that has mean τ and an arbitrarily small variance. Recall that the sum of k independent exponential

random variables with rate r is gamma distributed with rate parameter r, shape parameter k, mean

k/r, variance k/r2, and coefficient of variation 1/
√
k. Thus, fixing the mean at τ by setting the rate

r = k/τ, it follows that increasing k can yield an arbitrarily small coefficient of variation.

We may therefore use the Erlang case of Corollary 1.1 to write down a system of ODEs that

approximates eq. (1.5). By the above (and see [34] and references therein), it follows that eq. (1.5) can

be approximated by

d

dt
x0(t) = f(x0(t),xn(t)) (1.6a)

d

dt
xi(t) =

k

τ
[xi−1(t)−xi(t)], i = 1, . . . ,k (1.6b)

with initial conditions xi(0) = φ(−τ i/k) for i = 0, . . . ,k. Here, x0(t) ≈ x(t) and that approximation
improves as the value of k increases.

2. Results

With the above preliminaries in hand, we may now detail our GLCT-based procedure which serves as

a tool for implementing the GLCT to derive new ODE models with phase-type dwell time assumptions.

2.1. GLCT-based Procedure For Deriving ODE Models. Many existing mean field state transi-

tion models, in the form of an ODE, DDE, or distributed delay equation can be extended by a system

of ODEs derived using the following GLCT-based procedure.

(1) For an ODE model (with exponential dwell times), it is usually straightforward to obtain a

generalized model using Lemma 1.1 directly. In the case of a distributed delay equation or a

DDE model, this first step entails generating a system of ODEs using the standard linear chain

trick as detailed in section 1.2 above.

(2) The ODE model obtained in the first step above can then be written in matrix-form, as suggested

by the form of equations in Lemma 1.1 and Corollary 1.1. The resulting set of equations will

then be in a matrix form, where the matrix-vector pairs that were parameterized for the Erlang

distributions used to derive the model (see eq. (1.4)) can now be replaced by a more general

vector and matrix pair that together parameterize an arbitrary phase-type distribution.

In the sections below, we will illustrate the application of this GLCT-based procedure for deriving

new ODE models by generalizing various biological models taken from the peer reviewed literature. For

each model, we first introduce the application context and highlight some of the key model assumptions

related to delays and the distribution of times spent in different states. This entails viewing each

model as mean field equations corresponding to some unspecified stochastic model. We then derive

from each model an ODE model that incorporates phase-type distributed delays or dwell times, thereby

generalizing or approximating the original model.

2.2. Model 1: Tumor Growth Inhibition (TGI) Model. Simeoni et al [32] introduced a simple

model of tumor growth inhibition that employs the standard Linear Chain Trick (LCT) to incorporate

an Erlang distributed time to cell death following tumor cell damage from treatment. That model was

subsequently analyzed using standard approaches from dynamical systems [21, 20], and has been used

elsewhere in the study of tumor growth and the development of cancer treatments [29, 33]. The Simeoni

model has also been extended to a ‘‘competing Poisson processes” type of assumption (compare Fig. 6



416 P. HURTADO AND C. RICHARDS

Z0 Z1 Z2 Z3
cell 

death
k0 c(t) k1 k1 k1

Tumor
Growth

GF(z0,w)

Figure 1. Schematic diagram of the Tumor Growth Inhibition model (TGI) by Simeoni

et al [32]. See the main text for further details.

in [16] to Fig. 2 in [29] and Fig. 1 in [36]) in order to model tumor cell death arising from the combined

effects of two drugs with no pharmacokinetic interaction [29, 36].

The basic TGI model in [32] is given by eqs. (2.1) and (2.2) below. In the absence of pharmacological

treatment, the amount of cycling (replicating) tumor cells at time t, z0(t), grows according to the overall

growth rate2

GF(z0,w) =
λ0 z0(t)[

1 +

(
λ0
λ1
w(t)

)ψ] 1
ψ

. (2.1)

Treatment is assumed in [32] to begin at time t0 > 0, and accordingly the effect of that treatment

c(t) = 0 for 0 ≤ t ≤ t0. Once treatment begins, cells that are damaged by the treatment then progress
through a series of states Zi, i = 1, . . . ,n, prior to cell death (see Fig. 1). Together, the full model is

given in [32] by

dz0(t)

dt
= GF(z0(t),w(t))−k0 c(t)z0(t) (2.2a)

dz1(t)

dt
= k0 c(t)z0(t)−k1 z1(t) (2.2b)

dzi(t)

dt
= k1 zi−1(t)−k1 zi(t), i = 2, . . . ,n (2.2c)

w(t) =

n∑
i=0

zi; z0(0) = w0, zi(0) = 0, i = 1, . . . ,n (2.2d)

where w is the total amount of tumor cells, and k0 and c(t) ≥ 0 determine the rate of initial tumor
cell damage from the treatment.

The above assumptions and interpretations from [32] can also be described as follows. Viewed through

the lens of the GLCT, k0 and c(t) determine the distribution of time spent in the base state Z0, which

follows the first event time distribution under a non-homogeneous Poisson process with rate r(t) = k0c(t)

(for details, see[16]). Parameters n and k1 are the shape and rate parameters, respectively, for the Erlang

distributed time until cell death for the cells damaged by the treatment. The treatment is assumed to

have no effect on the time until cell death after the initial damage to the cell.

2.2.1. Generalized TGI Model. The GLCT-based procedure described in section 2.1 can be used to

extend this model to instead assume a more general phase-type distributed time to cell death. The

2This growth rate function is an approximation of the piece-wise function that is equal to λ0z0 when w < λ0/λ1, and

λ1z0/w when w ≥ λ0/λ1. See [21] for details.



A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 417

equations for zi, i = 1, . . . ,n in eqs. (2.2) can be written in matrix form, using Lemma 1.1, where

v =



1

0
...

0


 and M =




−k1 k1 0 · · · 0

0 −k1 k1
... 0

...
...

...
...

...

0 0
... −k1 k1

0 0 · · · 0 −k1



. (2.3)

This yields the more compact, and more general, set of equations below,

dz0(t)

dt
= GF(z0(t),w(t))−k0 c(t)z0(t) (2.4a)

dx(t)

dt
= k0 c(t)z0(t) v + M

Tx (2.4b)

where x = [z1,z2, . . . ,zn]
T, w(t) =

∑n
i=0 zi(t), z0(0) = w0, zi(0) = 0, i = 1, . . . ,n.

Note that eqs. (2.4) generalize the TGI model in the sense that these equations accommodate any

phase-type distribution assumption for the time to cell death following the initial effect of treatment,

not just the Erlang distribution assumed in the original TGI model and parameterized by eqs. (2.3).

Additionally, this matrix-vector form of the original TGI model (i.e., assuming an Erlang distribution)

can still be used with some benefit for both computational and mathematical analyses of the TGI model

given by eqs. (2.2), where those analyses can take advantage of the matrix-vector form of these more

general equations.

2.3. Model 2: Perscription Opioid Epidemic Model. Consider the prescription opioid epidemic

model by Battista et al [2], which is a system of ordinary differential equations with no explicit time

delays, and (implicit) exponentially distributed dwell times in multiple states.

The model assumes individuals are in one of four different states: S, P , A, and R. Here S respresents

the susceptible class. These individuals are not using opioids or recovering from addiction. P represents

the prescribed users (those who are prescribed the drugs and using them but have no addiction). A

represents addicted individuals who can be using either prescribed or ilicit opioids, and R represents the

class of individuals undergoing treatment and rehabilitation to recover from addiction.

The model as given in [2], where the dot over each state variable indicates a time derivative, is

Ṡ = −αS −βASA−βPSP + �P + δR + µ(P + R) + µ∗A (2.5a)

Ṗ = αS − (� + γ + µ)P (2.5b)

Ȧ = γP + σR + βASA + βPSP − (ζ + µ∗)A (2.5c)

Ṙ = ζA− (δ + σ + µ)R. (2.5d)

The term αS is the number of individuals transitioning from the susceptible state to the prescribed

state after being prescribed opioids per unit time, βASA is the rate of those transitioning from state S to

state A after interacting with addicted individuals, and similarly βPSP represents the rate of individuals

who transition from the susceptible class to the addicted class after exposure to opioids via perscription

opiod users who have extra or unsecured drugs. The terms �P and δR are the rate individuals leave the

prescribed users class without becoming addicted and then reenter the susceptible class at per-capita

rate �, and those who leave the rehabilitation state after treatment and reenter the susceptible state

at per-capita rate δ. The rates µP, µR and µ∗A are the death rates for the prescribed, rehabilitated,

and addicted classes (to ensure constant population size, deaths are replaced instantaneously by new

susceptible individuals). The term γP is the rate that individuals leave the prescribed class by becoming



418 P. HURTADO AND C. RICHARDS

addicted to their prescription opioids, ζA is the rate at which addicted individuals initiate treatment,

and σR is the rate at which individuals who are undergoing treatment reenter the addiction class.

We may also interpret the model terms described above as follows. Focusing on eq. (2.5b), for

example, prescription users remain in the prescribed state for an exponentially distributed amount time

(with rate � + γ + µ), and the proportion of individuals which leave the prescribed user state, and go

to the susceptible state, is �
�+γ+µ

(see section 3.5.5 in [16]). It follows that the net rate of individual

entering the susceptible state from the prescribed state is therefore �
�+γ+µ

(� + γ + µ)P = �P . Similarly,

the proportion of individuals who go on to become addicted, and who die, are given by γ
�+γ+µ

and
µ

�+γ+µ
, respectively.

2.3.1. Generalized Opioid Epidemic Model. We can use the procedure described in section 2.1 to extend

this model by replacing the implicit assumption of exponentially distributed dwell times in each state

with arbitrary phase-type distributed dwell times.

Assume the dwell time distribution for the prescribed user state P is a continuous phase-type

distribution parameterized by the n× 1 parameter vector vP and n×n matrix MP. Then to total
number of individuals in state P is given by the sum of the n sub-states Pi, i = 1, . . . ,n. Let

x = [P1,P2, . . . ,Pn]
T
. Then by the GLCT (Lemma 1.1), the mean field equations for our prescribed

user sub-states are

ẋ = vPαS + MP
Tx.

Observe that if we let vP be a one dimensional row vector with its first and only entry being a 1 and

let MP = [−(� + γ + µ)] be a 1×1 matrix, we will arrive at our original equation, eq. (2.5b).
Recall that individuals who leave the prescribed user state either transition to the addicted state, the

susceptible state, or they die. We can denote these proportions as FPA, FPS, and FPD, respectively, where

Fij ∈ [0,1] and FPA + FPS + FPD = 1. Note that in the original model FPA = γ(�+γ+µ),FPS =
�

(�+γ+µ)
,

and FPD =
µ

(�+γ+µ)
. This yields the model:

Ṡ = −αS −βASA−βPSP + FPS(−MP1)Tx + δR + FPD(−MP1)Tx + µR + µ∗A (2.6a)

ẋ = vPαS + MP
Tx (2.6b)

Ȧ = FPA(−MP1)Tx + σR + βASA + βPSP − (ζ + µ∗)A (2.6c)

Ṙ = ζA− (δ + σ + µ)R. (2.6d)

Similarly, we can generalize the addicted and rehabilitated states with phase-type dwell time distri-

butions, assuming the respective phase-type distributions are parameterized by vA, MA, vR, and MR.

Let y = [A1,A2, . . . ,Ak]
T
denote the k sub-states of A, and z = [R1,R2, . . . ,Rm]

T
the m sub-states of

R. This yields the generalized model:

Ṡ = −αS −βASA−βPSP + (FPS + FPD)(−MP1)Tx +

(FRS + FRD)(−MR1)T z + FAD(−MA1)T y
(2.7a)

ẋ = vPαS + MP
Tx (2.7b)

ẏ = vA
(
FPA(−MP1)Tx + FRA(−MR1)Tz + βASA + βPSP

)
+ MA

T y (2.7c)

ż = vR(FAR(−MA1)T y) + MRTz. (2.7d)

It is worth noting that the original model eqs. (2.5) from [2] are a special case of eqs. (2.7), as are

any intermediate extensions of the original model obtained by applying the standard Linear Chain Trick

(LCT) to impose Erlang distributed dwell times on one or more of the four main states.



A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 419

2.4. Model 3: Within-Host Model of The Immune-Pathogen Interaction. In [15], Hurtado

incorporated a specific (adaptive) immune response component to the innate immune response model by

Reynolds et al [28]. The scaled version of this within-host model, as stated in [15], is

dp

dt
= kpgp(1−p)−

kmp

µp + p
−K(y)np (2.8a)

dn

dt
=

n + kp p

xn + n + kp p
−µn p (2.8b)

dy0
dt

=
(np)α

xαy + (np)
α
−µy0 y0 (2.8c)

dy

dt
= µy0 y0 −µy y (2.8d)

In this model, p is the scaled pathogen (bacteria) population size, which follows a logistic growth model

in the absence of an immune response. The second term in eq. (2.8a) models the effect of some baseline

local immune defenses capable of neutralizing a small population of pathogen, and mathematically

introduces a strong Allee effect into the model. The level of innate immune activity n increases in

response to the presence of pathogen, as well as from a positive feedback loop, and the interaction

of this innate immune activity and pathogen stimulates progenitor cells (y0) that mature into active

specific immune components (y), e.g., B-cells, which augment the pathogen-killing capacity of the innate

immune components (i.e., which increase K(y)). For further details on this model, see [15] and [28].

In this model, the specific immune response delay can be thought of as an exponentially distributed

maturation time (with mean 1/µy0) and the duration of the active immune response (i.e., the dwell time

of mature specific immune components in state y) is also exponentially distributed (with mean 1/µy).

2.4.1. Within-Host Model With Phase-Type Delays In The Specific Immune Response. Both of these

exponential dwell time distribution assumptions associated with the specific immune response can be

replaced by phase-type distributions with respective parameters vy0, My0, vy and My, respectively. To

do this, we first partition state Y0 into sub-states Xi, i = 1, . . . ,m, and the state Y into sub-states Zj,

j = 1, . . . ,n, where y0 =
∑m
i=1 xi and y =

∑n
j=1 zj, and we let x = [x1, . . . ,xm]

T and z = [z1, . . . ,zn]
T.

The GLCT (Lemma 1.1) then yields the more general model

dp

dt
= kpgp(1−p)−

kmp

µp + p
−K(y)np (2.9a)

dn

dt
=

n + kp p

xn + n + kp p
−µn p (2.9b)

dx

dt
= vy0

(np)α

xαy + (np)
α
+ My0

T x (2.9c)

dz

dt
= −1TMyTvy + MyTz. (2.9d)

2.5. Model 4: Cell-To-Cell Spread of HIV. In [5], Culshaw et al introduce an integro-differential

equation model of the cell-to-cell spread of HIV, which incorporates a distributed time delay in the

time between cells becoming infected and infectious. They then derive from this general model multiple

other models which differ only in the specific assumptions on the form of this delay distribution. The

consider models with no delay, with a fixed (Dirac-δ) delay of τ time units, and with an exponentially

distributed delay. Here we will extend these models to a general phase-type distributed delay.

In the most general model described in [5], state variable C(t) represents the concentration of healthy

cells at time t, and I(t) is the concentration of infected cells, and



420 P. HURTADO AND C. RICHARDS

dC

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.10a)

dI

dt
= k′i

∫ t
−∞

C(u)I(u)F(t−u)du−µII(t). (2.10b)

Parameter rC is the net growth rate of the healthy cell population, CM is an effective carrying

capacity of the system, kI is an infection rate parameter, k
′
I/kI is the fraction of cells surviving the

incubation period, and µI is the per capita death rate of infected cells (implicitly, the infected cell

lifetime is exponentially distributed with mean 1/µI). Initial values for C and I must be functions

defined over all s ∈ (−∞,0] and are denoted φ(s) ≥ 0 and ψ(s) ≥ 0, respectively.
In [5], the delay kernel F(u) is assumed to be a step function from 0 up to 1 at u = τ ≥ 0, or (in the

exponential case) of the form

F(u) =
αn+1un

n!
e−αu (2.11)

which is just the density function for an Erlang distribution with rate α and shape n + 1 (and thus,

mean (n + 1)/α and coefficient of variation 1/
√
n + 1). The weak and strong kernels referenced in [5]

are just the particular cases where the shape parameter is 1 (i.e., an exponential distribution with rate

α) or 2 (Erlang with rate α and shape 2), respectively.

In [5], Culshaw et al derive the specific models described above from this more general integro-

differential equation model (eqs. (2.10)), which we have summarized below (see [5] for a comparison of

the dynamics of these three models).

First, assuming F(u) = δ(u − τ) is a Dirac-δ function at time τ ≥ 0 yields the delay differential
equation (or ODE if τ = 0) below, as written in [5], with the same initial conditions as eqs. (2.10) if

τ > 0 or with initial conditions C(0) = c0 ≥ 0 and I(0) = I0 ≥ 0 if τ = 0.

dC(t)

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.12a)

dI(t)

dt
= k′iC(t− τ)I(t− τ)−µII(t). (2.12b)

Next, the authors in [5] assumed a ‘‘weak kernel” (i.e., exponentially distributed delay with rate α)

and derived the following system of ODEs:

dC(t)

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.13a)

dX(t)

dt
= αC(t)I(t)−αX(t) (2.13b)

dI(t)

dt
= k′I X(t)−µII(t). (2.13c)

2.5.1. ODE Model of Cell-To-Cell HIV Spread With Phase-Type Lags In Cells Becoming Infectious.

We can extend the above models as follows. First, substituting Y (t) = kI
α
X(t) yields a more natural (in



A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 421

the context of the LCT, and in terms of the units of X and Y ) set of equations

dC(t)

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.14a)

dY (t)

dt
= kI C(t)I(t)−αY (t) (2.14b)

dI(t)

dt
=
k′I
kI
αY (t)−µII(t). (2.14c)

From these equations, it is straightforward to derive a more general ODE model using our GLCT-based

procedure. Applying the LCT to eqs. (2.14) yields eqs. (2.15) below, which correspond to any choice

of α > 0 and non-negative integer value n ≥ 0 for the delay kernel F in eq. (2.11) (i.e., any Erlang
distribution with shape n + 1 and rate α).

dC(t)

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.15a)

dY1(t)

dt
= kI C(t)I(t)−αY1(t) (2.15b)

dYi(t)

dt
= αYi−1(t)−αYi(t), i = 2, . . . ,n + 1 (2.15c)

dI(t)

dt
=
k′I
kI
αYn+1(t)−µII(t). (2.15d)

Re-writing the above equations in the particular matrix form suggested by the GLCT (Lemma 1.1)

yields the more general set of equations below, which are the desired set of model equations for which

the Erlang distribution assumption (with parameters n + 1 and α) have been replaced by a phase-type

distribution parameterized by the length k vector v and k ×k matrix M, where y = [Y1, . . . ,Yk]T.
dC(t)

dt
= rCC(t)

(
1−

C(t) + I(t)

CM

)
−kIC(t)I(t) (2.16a)

dy(t)

dt
= kI C(t)I(t) v + M

T y(t) (2.16b)

dI(t)

dt
= −

k′I
kI

1T MT y(t)−µII(t). (2.16c)

Further generalizations, e.g., to accommodate time-varying v, M, or survival fraction f = k′I/kI are

also possible, as detailed in [16]. Note that it is also straightforward to further extend this model by

using Lemma 1.1 to replace the exponential dwell time assumption (with rate µI) in state I, eq. (2.16c)

by a phase-type distribution assumption.

3. Discussion

Mean field state transition models, written as ordinary differential equations (ODEs), are widely

used throughout the sciences, but too often they include overly simplistic assumptions regarding time

delays and the duration of time individual entities spend in specific states. In this paper, we introduce

a relatively simple, novel procedure for implementing the generalized linear chain trick (GLCT) to

expedite the derivation of ODE models that can incorporate a broader range of underlying assumptions.

We use this procedure to derive new models, based upon existing models taken from the peer reviewed

literature, illustrating the relative ease of deriving new models within the context of the GLCT. These

examples include a mix of biological applications and different mathematical frameworks (DDEs, ODEs,

and integro-differential equations), which reflect a mix of different implicit and explicit delay and dwell

time assumptions.



422 P. HURTADO AND C. RICHARDS

These straightforward generalizations illustrate how modelers can easily incorporate phase-type

distributed delays and dwell times into ODE models, and how underlying stochastic model assumptions

are reflected in the structure of corresponding ODE models. Importantly, these more general model

formulations can also be used in the computational and mathematical analysis of models that only

assume Erlang distributions (i.e., that could otherwise be derived using the standard linear chain trick).

For an example, in [17] we illustrate the potential computational benefits of using a GLCT formulation

of models with Erlang dwell time assumptions when computing numerical solutions to such models.

The benefits of using ODE models derived using the GLCT extend beyond the benefits of writing the

model in matrix form. For example, concepts from Markov chain theory, including reward process

theory, can be employed in the analysis (and interpretation of results) of ODE models derived using

this GLCT-based procedure [18].

These generalized models also lay the groundwork for data-driven model formulations that incorporate

Coxian, hyperexponential, hypoexponential (i.e, generalized Erlang) and other phase-type distributions

into these and similar mean field ODE models. Statistical tools such as BuTools [14, 13] allow modelers

to fit phase-type distributions to data, thereby allowing modelers to build approximate empirical

distributions into ODE models using the GLCT. However, it is important to note that there are

limitations to approximating some delay or dwell time assumptions with phase-type distributions. For

example, delay distributions with compact support (e.g., a continuous uniform distribution) may not

be well approximated by phase-type distributions. There are also models with more complex state

transition or dwell time assumptions that may be more appopriate to derive using the more rigorous

approaches detailed in [16].

In closing, we hope this GLCT-based procedure for the derivation of new ODE models – that

approximate or generalize existing models – proves to be helpful to modelers in their efforts to build

better models, to check the consequences of certain simplifying assumptions, and to gain better intuition

for how underlying assumptions are reflected in the structure of ODE model equations.

Acknowledgements

The authors thank the anonymous reviewers of this manuscript for their constructive feedback, and Dr.

Deena Schmidt for conversations and comments that improved this manuscript. The authors also thank

the organizers and sponsors of the Second International Conference on Applications of Mathematics

to Nonlinear Sciences (AMNS-2019), held June 27-30, 2019, in Pokhara, Nepal, and the editors of this

Thematic Issue.

Funding

This work was supported by a grant awarded to PJH by the Sloan Scholars Mentoring Network

of the Social Science Research Council with funds provided by the Alfred P. Sloan Foundation; and

this material is based upon work supported by the National Science Foundation under Grant No.

DEB-1929522.

Disclosure statement

The authors declare that they have no conflict of interest.

References

1. R.M. Anderson and R.M. May, Infectious Diseases of Humans: Dynamics and Control, Oxford University Press,

1992.

2. N.A. Battista, L.B. Pearcy, and W.C. Strickland, Modeling the prescription opioid epidemic, B. Math. Biol. 81 (2019),

no. 7, 2258--2289.



A PROCEDURE TO DERIVE NEW ODE MODELS VIA THE GLCT 423

3. A. Beuter, L. Glass, M.C. Mackey, and M.S. Titcombe (eds.), Nonlinear Dynamics in Physiology and Medicine,

Interdisciplinary Applied Mathematics (Book 25), Springer, September 2003.

4. M. Bladt and B.F. Nielsen, Phase-type distributions, Matrix-Exponential Distributions in Applied Probability, Springer

US, 2017, pp. 125--197.

5. R.V. Culshaw, S. Ruan, and G. Webb, A mathematical model of cell-to-cell spread of HIV-1 that includes a time

delay, J. Math. Biol. 46 (2003), no. 5, 425--444.

6. O. Diekmann, M. Gyllenberg, and J.A.J. Metz, Finite dimensional state representation of linear and nonlinear delay

systems, J. Dyn. Differ. Equ. (2017), 1439--1467.

7. O. Diekmann and J.A.P. Heesterbeek, Mathematical epidemiology of infectious diseases: Model building, analysis and

interpretation, Wiley Series in Mathematical and Computational Biology, John Wiley & Sons, LTD, New York, 2000.

8. Z. Feng, D. Xu, and H. Zhao, Epidemiological models with non-exponentially distributed disease stages and applications

to disease control, B. Math. Biol. 69 (2007), no. 5, 1511--1536.

9. Z. Feng, Y. Zheng, N. Hernandez-Ceron, H. Zhao, J.W. Glasser, and A.N. Hill, Mathematical models of Ebola---

Consequences of underlying assumptions, Math. Biosci. 277 (2016), 89--107.

10. T. Gedeon and G. Hines, Upper semicontinuity of morse sets of a discretization of a delay-differential equation, J.

Differ. Equations 151 (1999), no. 1, 36--78.

11. , Upper semicontinuity of morse sets of a discretization of a delay-differential equation: An improvement, J.

Differ. Equations 179 (2002), no. 2, 369--383.

12. A. Horváth, M. Scarpa, and M. Telek, Phase type and matrix exponential distributions in stochastic modeling,

Principles of Performance and Reliability Modeling and Evaluation: Essays in Honor of Kishor Trivedi on his 70th

Birthday (Lance Fiondella and Antonio Puliafito, eds.), Springer International Publishing, Cham, 2016, pp. 3--25.

13. G. Horváth and M. Telek, Butools v2.0, http://webspn.hit.bme.hu/~telek/tools/butools/doc/, Accessed: May 15,

2020.

14. , BuTools 2: A rich toolbox for Markovian performance evaluation, ValueTools 2016 - 10th EAI International

Conference on Performance Evaluation Methodologies and Tools, Association for Computing Machinery, 1 2017,

pp. 137--142 (English).

15. P.J. Hurtado, Within-host dynamics of mycoplasma infections: Conjunctivitis in wild passerine birds, J. Theor. Biol.

306 (2012), 73 -- 92.

16. P.J. Hurtado and A.S. Kirosingh, Generalizations of the ‘linear chain trick’: incorporating more flexible dwell time

distributions into mean field ODE models, J. Math. Biol. 79 (2019), no. 5, 1831--1883.

17. P.J. Hurtado and C. Richards, Building mean field state transition models using the generalized linear chain trick

and continuous time markov chain theory, arXiv:2007.03902, 2020.

18. , Finding reproduction numbers for epidemic models & predator–prey models of arbitrary finite dimension

using the generalized linear chain trick, arXiv:2008.06768, 2020.

19. C. Lin, L. Wang, and G.S.K. Wolkowicz, An alternative formulation for a distributed delayed logistic equation, B.

Math. Biol. 80 (2018), no. 7, 1713--1735.

20. P. Magni, M. Germani, G. De Nicolao, G. Bianchini, M. Simeoni, I. Poggesi, and M. Rocchetti, A minimal model of

tumor growth inhibition, IEEE. Trans. Biomed. Eng. 55 (2008), no. 12, 2683--2690.

21. P. Magni, M. Simeoni, I. Poggesi, M. Rocchetti, and G. De Nicolao, A mathematical model to study the effects of

drugs administration on tumor growth dynamics, Math. Biosci. 200 (2006), no. 2, 127--151.

22. J.A.J. Metz and O. Diekmann (eds.), The dynamics of physiologically structured populations, Lecture Notes in

Biomathematics, vol. 68, Springer, Berlin, Heidelberg, 1986.

23. , Exact finite dimensional representations of models for physiologically structured populations. I: The abstract

formulation of linear chain trickery, Proceedings of Differential Equations With Applications in Biology, Physics, and

Engineering 1989 (J. A. Goldstein, F. Kappel, and W. Schappacher, eds.), vol. 133, 01 1991, pp. 269--289.

24. W.T. Mocek, R. Rudnicki, and E.O. Voit, Approximation of delays in biochemical systems, Math. Biosci. 198 (2005),

no. 2, 190--216.

25. W.W. Murdoch, C.J. Briggs, and R.M. Nisbet, Consumer--resource dynamics, Monographs in Population Biology,

vol. 36, Princeton University Press, Princeton, USA, 2003.

26. R.M. Nisbet, W.S.C. Gurney, and J.A.J. Metz, Stage structure models applied in evolutionary ecology, Applied

Mathematical Ecology (Simon A. Levin, Thomas G. Hallam, and Louis J. Gross, eds.), Springer Berlin Heidelberg,

Berlin, Heidelberg, 1989, pp. 428--449.

27. P. Reinecke, L. Bodrog, and A. Danilkina, Phase-Type Distributions, Resilience Assessment and Evaluation of

Computing Systems (Katinka Wolter, Alberto Avritzer, Marco Vieira, and Aad van Moorsel, eds.), Springer Berlin

Heidelberg, Berlin, Heidelberg, 2012, pp. 85--113.

http://webspn.hit.bme.hu/~telek/tools/butools/doc/


424 P. HURTADO AND C. RICHARDS

28. A. Reynolds, J. Rubin, G. Clermont, J. Day, Y. Vodovotz, and G.B. Ermentrout, A reduced mathematical model

of the acute inflammatory response: I. Derivation of model and analysis of anti-inflammation, J. Theor. Biol. 242

(2006), no. 1, 220 -- 236.

29. M. Rocchetti, F. Del Bene, M. Germani, F. Fiorentini, I. Poggesi, E. Pesenti, P. Magni, and G. De Nicolao, Testing

additivity of anticancer agents in pre-clinical studies: A PK/PD modelling approach, Eur. J. Cancer 45 (2009), no. 18,

3336--3346.

30. M.R. Roussel, The use of delay differential equations in chemical kinetics, J. Phys. Chem. 100 (1996), no. 20,

8323--8330.

31. S. Ruan, Delay differential equations in single species dynamics, Delay Differential Equations and Applications

(O. Arino, M.L. Hbid, and E. Ait Dads, eds.), NATO Science Series, vol. 205, Springer, 2006, pp. 477--517.

32. M. Simeoni, P. Magni, C. Cammia, G. De Nicolao, V. Croci, E. Pesenti, M. Germani, I. Poggesi, and M. Roc-

chetti, Predictive pharmacokinetic-pharmacodynamic modeling of tumor growth kinetics in xenograft models after

administration of anticancer agents, Cancer Res. 64 (2004), no. 3, 1094--1101.

33. M. Simeoni, G. De Nicolao, P. Magni, M. Rocchetti, and I. Poggesi, Modeling of human tumor xenografts and dose

rationale in oncology, Drug Discovery Today: Technologies 10 (2013), no. 3, e365--e372.

34. H. Smith, An introduction to delay differential equations with applications to the life sciences, vol. 57, Springer

Science & Business Media, 2010.

35. S.H. Strogatz, Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering, 2nd

ed., Studies in Nonlinearity, Westview Press, July 2014.

36. N. Terranova, M. Germani, F. Del Bene, and P. Magni, A predictive pharmacokinetic–pharmacodynamic model of

tumor growth kinetics in xenograft mice after administration of anticancer agents given in combination, Cancer

Chemoth. Pharm. 72 (2013), no. 2, 471--482.

37. H.J. Wearing, P. Rohani, and M.J. Keeling, Appropriate models for the management of infectious diseases, PLOS

Med. 2 (2005), no. 7, 0621--0627.

Corresponding Author, Department of Mathematics and Statistics, University of Nevada, Reno; Reno,

Nevada, USA

E-mail address: phurtado@unr.edu

Department of Mathematics and Statistics, University of Nevada, Reno; Reno, Nevada, USA

E-mail address: cj.richards@nevada.unr.edu


	1. Introduction
	1.1. Generalized Linear Chain Trick
	1.2. Deriving ODEs from Delay Differential Equations & Distributed Delay Equations

	2. Results
	2.1. GLCT-based Procedure For Deriving ODE Models
	2.2. Model 1: Tumor Growth Inhibition (TGI) Model
	2.3. Model 2: Perscription Opioid Epidemic Model
	2.4. Model 3: Within-Host Model of The Immune-Pathogen Interaction
	2.5. Model 4: Cell-To-Cell Spread of HIV

	3. Discussion
	Acknowledgements
	Funding
	Disclosure statement
	References