\begin{matrix}c \\ ◯ \end{matrix} Revista de Matemática: Teoría y Aplicaciones 2020 27(1) : 49–71 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online)

doi: https:// doi. org/ 10.15517/ rmta. v27i1. 39948

a delay differential equations model for disease transmission dynamics

un modelo de ecuaciones diferenciales con retraso para la dinámica de transmisión de enfermedades

Mustafa Erdem^{∗}  Muntaser Safan^{†}

Carlos Castillo-Chavez^{‡}

Received: 18/May/2019; Revised: 20/6/2019; Accepted: 17/Sep/2019

∗Arizona State University, Simon A. Levin Mathematical, Computational and Modeling Sciences Center, Tempe AZ, United States. E-Mail: mustafa.erdem@asu.edu

†Arizona State University, Simon A. Levin Mathematical, Computational and Modeling Sciences Center, Tempe, United States. Mansoura University, Mathematics Department, Faculty of Science, Mansoura, Egypt. Umm Al–Qura University, Department of Mathematical Sciences, Faculty of Applied Sciences, Makkah, Saudi Arabia. E-Mail: muntaser_safan@yahoo.com

‡Arizona State University, Simon A. Levin Mathematical, Computational and Modeling Sciences Center, Tempe AZ, United States; Brown University, Visiting Provost Professor of Applied Mathematics. E-Mail: carlos_castillo-chavez@brown.edu

Abstract

A delay differential equations epidemic model of SIQR (Susceptible-Infective-Quarantined-Recovered) type, with arbitrarily distributed periods in the isolation or quarantine class, is proposed. Its essential mathematical features are analyzed. In addition, conditions that support the existence of periodic solutions via Hopf bifurcation are identified. Nonexponential waiting times in the quarantine/isolation class lead not only to oscillations but can also support stability switches.

Keywords: delay differential equation; integro-differential equation; epidemic model; quarantine; stability switch; oscillations; stage structure.

Resumen

Se propone un modelo epidémico de ecuaciones diferenciales con retraso del tipo SIQR (por sus siglas en inglés) (Susceptible-Infeccioso-En cuarentena-Recuperado), con períodos arbitrariamente distribuidos en la clase de aislamiento o cuarentena. Se analizan sus características matemáticas esenciales. Además, se identifican las condiciones que respaldan la existencia de soluciones periódicas a través de la bifurcación de Hopf. Los tiempos de espera no exponenciales en la clase de cuarentena/aislamiento conducen no solo a oscilaciones sino que también pueden soportar cambios de estabilidad.

Palabras clave: ecuación diferencial con retraso; ecuación integro-diferencial; modelo epidémico; cuarentena; cambio de estabilidad; oscilaciones; estructura por etapas.

Mathematics Subject Classification: 92D25, 92D30, 92B99.

1 Introduction

The emergence of SARS in 2003 and the global impact of the recent 2009 influenza A (H1N1) pandemic have re-invigorated (for a historical perspective see[16]) our interests in the systematic use of policies of isolation, quarantine (ı & q), and/or social distancing as population-level methods of disease control [7, 26, 18]. The assessment, from a public health policy perspective, of the effectiveness of ı & q intervention measures is tied into their ability to reduce transmission at the population-level. A relevant set of questions include: How many individuals must be quarantined or isolated to significantly alter the transmission disease patterns at the population level and at what cost? Some forms of ı & q interventions have proved effective in the case of SARS [6].

However, from an economic standpoint, the implementation of ı & q policies is not cost-effective for all infectious diseases [22, 26, 18, 14]. The use of ı & q in the case of rubella in China, for example, may lead to higher disease levels since under the current system of no control, about 97% of the population acquires rubella antibodies from direct exposure to infectious individuals. China’s high levels of naturally-induced coverage would be impossible to match if, for example, China chose to put in place the highly effective US or Canada vaccination policies (see [15]).

Mathematical modeling and their mathematical analysis help assess the likelihood that ı & q are effective disease management strategies. In fact, several modeling studies have been carried out to study their effectiveness under the assumption that the waiting times in the ı or q classes are exponentially distributed, a most likely adequate assumption [22, 7]. Theoretical studies are routinely used to test the robustness of modeling assumptions. In this manuscript, we proceed to evaluate the robustness of a modeling critical assumption that is common to disease management models, namely that waiting times in the ı or q classes (here used interchangeably) are non-exponentially distributed. The results of the mathematical analysis are revealing. It is shown that the use of non-exponential waiting times in the quarantine/isolation class leads not only to oscillations but can also support stability switches. Whether or not this conclusion can be explicitly assessed in the context of specific naturally occurring epidemic outbreaks would depend on the availability of data. However, what the results in this manuscript show is that the massive implementation of such policies (as may be the case in pandemic situations or in deliberately engineered epidemics) could significantly impact the character of the observed disease incidence patterns.

The systematic mathematical description of the role of quarantine/isolation and the assessment of these interventions on disease dynamics has been studied [23, 5, 4, 11, 20, 9]. The models and results in this manuscript are tied into previous work on the role of quarantine on disease dynamics [23, 11, 20, 28, 9]. The role of variable latent periods in models for TB has been studied in [10]. The authors show that TB models including arbitrarily distributed latent stage and those including exponentially distributed period of latency lead to similar qualitative behaviours. Mathematical models with arbitrarily distributed periods of infections have been studied by Feng et al. [13, 12].

We introduce a compartmental model consisting of four classes, susceptibles S, infected/infectious I, isolated Q and recovered R individuals. Susceptible individuals who acquire the disease move to the I-class and a fraction of I-individuals are isolated/quarantined, that is, they are moved to the Q class while the rest of I individuals transit into the R-class after recovery.

The disease is assumed to be essentially non-fatal and therefore disease-induced deaths are considered to be “minimal” and therefore ignored. Hence, it is assumed that the total population N = S + I + Q + R is roughly constant. Explicitly, it is assumed that the birth/death per capita rates are equal to μ and that all newborns are susceptible. The per-capita rates at which individuals move from I to Q and from I or Q to R classes are respectively, θ and γ under exponentiallydistributed waiting time s. In general, Q-individuals move to the R-class at variable rates (non-exponentially distributed waiting times). In this manuscript, the mathematical discussion is focused on two extreme cases for comparison purposes: the situation when individuals are retained for a fixed period of time and the exponentially-distributed waiting times case.

The incidence rate is modeled using a quarantine-adjusted incidence rate βSI/(N − Q) instead of standard incidence rate, βSI/N–a mechanism that is known to be capable of sustaining periodic solutions via a Hopf–bifurcation [11, 20, 28].

Figure 1: SIQR model for influenza.

2 Formulation of the model

In this section, we introduce a general S − I − Q − R epidemic model that incorporates arbitrary waiting times in its epidemiological classes. The model derivation is standard (see [3]), but worth outlining in this general setting. If it is assumed that f(s) denotes the number of people who remain at a certain stage, say G, at time s and if ψ denotes the per-capita removal rate from stage G, then the rate of change of f(s) can be described by

\frac{df(s)}{ds} = −ψf(s),  0 ≤ s < ∞,  f(0) = f_{0}.

Therefore,

\frac{f(s)}{f_{0}} = e^{−ψs}, for s ≥ 0,

or, in other words, e^{−ψs} (usually known as the survival function) can be interpreted as denoting the proportion of individuals who were at stage G at time s = 0 and are still at stage G at time s = s. Consequently,

F(s) = \{ \begin{matrix}1 − e^{−ψs} & if s ≥ 0, \\ 0 & if s < 0, \end{matrix}

can be interpreted as the “probability” of or the proportion of individuals leaving stage G over the time interval [0, s). The individuals in the population of interest are classified according to their epidemiological status: S(t) denotes the susceptible class at time t; I(t) the infected/infectious class at time t; Q(t) the quarantined/isolated class at time t; and R(t) the recovered class with N(t) = S(t)+I(t)+Q(t)+R(t). Individuals are transferred from class to class according to the transfer flow box provided in Figure 1. The waiting times tied to the I and Q compartments are modeled by the function P : [0, ∞) → [0, 1] wherein P (h) denotes the proportion of individuals that (if alive) remain h units of time in the infective class after acquiring an infection at time t = 0. That is, P (h) is the sojourn function, cf. [29, 27, 28]; P (h) is nonnegative, nonincreasing and P (0) = 1. In order to clarify the role of each waiting times distribution, we introduced the following expanded notation. We let P (μ, h), P (γ, h) and P (θ, h) denote: the fraction of individuals at time t = 0 who are still alive at time t = h; the fraction of individuals infected at time t = 0 that if alive have not recovered at time t = h; and the fraction of individuals infected at time t = 0 that have yet to be quarantined at time t = h, respectively. Under this notation, we have that P (h) = P (μ, h)P (γ, h)P (θ, h), where μ, γ and θ denote per-capita death rate, the per-capita recovery rate and per-capita quarantine rate respectively. Analogously, let L(μ, h) and L(α, h) denote: the fraction of quarantined individuals at time t = 0 who are still alive at time t = h and the fraction of individuals quarantined at time t = 0 that if alive have not recovered at time t = h, respectively. Then, L(h) = L(μ, h)L(α, h) denotes the proportion of individuals that, if alive, remain in the isolated class at time t = h, given that they were in this class at time t = 0.

Finally, if

B(t) := βS(t)\frac{I(t)}{N − Q(t)},

denotes the incidence rate (new cases of infection per unit time), then the following general nonlinear system of integro-differential equations is used to capture the dynamics of a communicable disease (like influenza, measles or common cold) in the presence of arbitrary waiting times in the I and Q classes.

S′(t) = μN − βS \frac{I}{N − Q} − μS,
I(t) = ∫_{−∞}^{t} B (x) P (μ, t − x) P (γ, t − x) P (θ, t − x) dx,  (1)
Q(t) = ∫_{−∞}^{t} ∫_{−∞}^{τ} B (x) P (μ, τ − x) P (γ, τ − x) (−P′ (θ, τ − x)) dxdτ

− ∫_{−∞}^{t}[∫_{−∞}^{s} ∫_{−∞}^{τ} B (x) P (μ, τ − x)P (γ, τ − x)(−P′ (θ, τ − x))

dxdτ](−L′(t − s)) ds,

where the′ means differentiation with respect to time t. Since 1−P (θ, t) denotes the proportion of individuals who leave the isolation class, then −P′(θ, t) and −L′(t) represent the rate at which individuals leave the infected class (moving to the isolation class) and the rate at which individuals leave the isolation class, by recovery or death, respectively. We will consider two types of waiting times in the Q class, exponential time of isolation and fixed time of isolation respectively in order to highlight two extremes.

2.1 Exponentially distributed isolation time

We assume that P (κ, t) and L(κ, t) are of the form e^{−κt}, that is, the cases of exponential times of survival (κ = μ), infection (κ = γ + θ) and isolation (κ = α) are considered. Under exponential waiting times, system (1) reduces to the following system of differential equations:

S′ = μN − βS \frac{I}{N − Q} − μS,
I′ = βS \frac{I}{N − Q} − (γ + θ + μ)I,  (2)
Q′ = θI − (α + μ)Q.

A detailed account of the dynamics of a system slightly more general than (2) can be found in [11, 20] and [9]. Feng and Thieme [11] studied an SIQR model in which everybody in infective class moved directly to the isolation class and then to the recovered class. They proved the existence of periodic solution via Hopf Bifurcation and estimated the period of the solutions using singular perturbation theory. Hethcote et al. [20] studied exactly the same model as (2) proving sustained oscillations (via Hopf bifurcation).

2.2 Fixed time of isolation

We now assume that individuals stay in the Q class a fixed period of time. That is, we model the waiting times in the quarantine class (Q) by

G(h) = \{ \begin{matrix}1 & if 0 ≤ h < ρ, \\ 0 & if ρ ≤ h, \end{matrix}

which corresponds to a fixed isolation period [3, 27]. The function L in (1) takes the form

L(h) = G(h)e^{−μh}.

If we replace the exponential waiting in the Q-class with a fixed period of time ρ, then −L′(h) represents the rate at which the individuals leave the class Q by recovery or death. The second term in the equation of Q(t) in (1) can be rewritten as

∫_{−∞}^{t} A (s) e^{−μ(t−s)} [μH (s − t + ρ) + δ (s − t + ρ)] ds,  (3)

where

A(s) = ∫_{−∞}^{s} ∫_{−∞}^{τ} B (x) P (μ, τ − x) P (γ, τ − x) [−P′ (θ, τ − x)] dxdτ,

H denotes the Heaviside function and δ denotes the Dirac Delta function. The evaluation of the integral in (3) gives

A(t − ρ)e^{−μρ} + μ ∫_{t−ρ}^{t} A(s)e^{−μ(t−s)} ds.  (4)

Combining (4) with the equation for Q(t) in (1) leads to the following expression for Q(t)

Q(t) = A(t) − A(t − ρ)e^{−μρ} − μ ∫_{t−ρ}^{t} A(s)e^{−μ(t−s)} ds,  (5)

or, equivalently,

Q(t) = θ ∫_{0}^{ρ} I(t − s)e^{−μs} ds.  (6)

Thus, the corresponding system of integro-differential equations is given by

S′(t) = μN − βS \frac{I}{N − Q} − μS,
I(t) = ∫_{−∞}^{t} B(x)e^{−(μ+γ+θ)(t−x)} dx,  (7)
Q(t) = θ ∫_{0}^{ρ} I(t − x)e^{−μx} dx,

or, equivalently, from system (7) we can derive the following system of delaydifferential equations

S′(t) = μN − βS(t)\frac{I(t)}{N − Q(t)} − μS(t),
I′(t) = βS(t)\frac{I(t)}{N − Q(t)} − (γ + θ + μ)I(t),  (8)
Q′(t) = θ (I(t) − I(t − ρ)e^{−μρ}) − μQ(t),

with S(0) = S_{0} > 0, Q(0) = Q_{0} > 0 and I(t) ≥ 0 on t ∈ [−ρ, 0].

The disease-free equilibrium of (8) is E_{0} = (N, 0, 0). Using linearization around E_{0}, it is not hard to see that this trivial equilibrium is locally asymptotically stable if and only if β < γ + θ + μ, that is when R_{0} < 1 and unstable when R_{0} > 1, where

R_{0} = \frac{β}{γ + θ + μ},  (9)

is the basic reproduction number. It is noteworthy that the basic reproductive number is independent of the residence time ρ in the class Q, as it should be.

The nontrivial equilibrium, E^{∗} = (S^{∗}, I^{∗}, Q^{∗}), which corresponds to the endemic steady state of System (8) is given by

S^{∗} = \frac{1 − Θ̂(1 − e^{−μρ})}{R_{0} − Θ̂(1 − e^{−μρ})}N,
I^{∗} = \frac{μ}{γ + θ + μ} · \frac{R_{0} − 1}{R_{0} − Θ̂(1 − e^{−μρ})}N,  (10)
Q^{∗} = Θ̂(1 − e^{−μρ}) · \frac{R_{0} − 1}{R_{0} − Θ̂(1 − e^{−μρ})}N,

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 49–71, Jan–Jun 2020

where

Θ̂ = \frac{θ}{γ + θ + μ}.

Remarkably, in the absence of isolation (i.e., θ = 0), the solutions (10) reduce to the nontrivial equilibrium components of the classic endemic SIR model [19]. Moreover, it is noticeable that the components (10) of the endemic equilibrium E^{∗} make biological sense if and only if the basic reproduction number R_{0} > 1. In other words, when R_{0} < 1, the equilibrium E^{∗} disappears.

3 Basic properties of the delay and integro-differential equation models

In this section, relevant mathematical results associated with systems (7) and (8) are obtained when individuals in the Q class are kept for a fixed period of time. In fact, we have established the following result:

Proposition 1 A solution of the system of integro-differential equations (7) satisfies the system of delay differential equations (DDE’s) (8). Conversely, if we let S(t), I(t), Q(t) be a solution of (8) with initial conditions given (on the intervals) as specified above and if in addition,

Q(0) = θ ∫_{−ρ}^{0} I(x)e^{μx} dx,  (11)

I(0) = ∫_{−∞}^{0} B(x)e^{kx} dx,  (12)

where k = θ + μ + γ, then this solution satisfies system (7). Moreover, for all t ≥ 0, the solution exists, is unique and has S(t) ≥ 0, I(t) ≥ 0 and Q(t) ≥ 0.

Proof. The assertion that the solution of the system of integro-differential equations (7) satisfies the system of the delay differential equations (DDE’s) (8) can be checked easily. Conversely, from (8) the equation for Q(t) it follows that

Q′(t) + μQ(t) = θI(t) − θI(t − ρ)e^{−μρ}.  (13)

Multiplying both sides by e^{μt}, we get

e^{μt} [Q′(t) + μQ(t)] = e^{μt} [θI(t) − θI(t − ρ)e^{−μρ}] ,  (14)

which on integrating leads to

e^{μt}Q(t) − Q(0) = ∫_{0}^{t} θe^{μτ} I(τ) dτ − ∫_{0}^{t} θI(τ − ρ)e^{μ(τ−ρ)} dτ.  (15)

Using change of variable x = τ − ρ results in

e^{μt}Q(t) = Q(0) − ∫_{−ρ}^{0} θI(x)e^{μx} dx + ∫_{t−ρ}^{t} θI(x)e^{μx} dx.  (16)

From (11), the equation for Q(t) in (7) follows. Similarly, by integrating the equation for I(t) in (8) and using (12), one arrives at the equation of I(t) in (7). The proof for the existence and uniqueness of solutions for System (8) can be established using standard theorems (see [17]), provided that the denominator N − Q(t) does not go to zero as t changes. The inequalities I(t) ≥ 0 and Q(t) ≥ 0 for t ≥ 0 follow immediately from the equations in (7) and (8) and the observation that, under the assumptions just made, I(t) ≥ 0 for t ∈ [−ρ, 0]. Consequently, the solution exists globally for all t > 0 and is unique.

Let us define the region D as:

D = {(S, I, Q) : S ≥ 0, I ≥ 0, Q ≥ 0, S + I + Q ≤ N} .

Proposition 2 If R_{0} ≤ 1, the disease free equilibrium, E_{0}, is globally asymptotically stable in the region D. Whenever R_{0} > 1, E_{0} is unstable and there exists a unique positive endemic equilibrium E_{1} given by (10).

Proof. The proof is carried out along the same line of the global stability proofs made in [20]. The Liapunov function V = I works out well here since

V′ = β \frac{SI}{N − Q} − (γ + θ + μ)I ≤ 0.

Therefore, by the Liapunov–LaSalle theorem and a straightforward application of [20, Lemma 1] we conclude that all the solutions in D should approach E_{0} whenever R_{0} ≤ 1.

In the following, we state another proposition that gives us an upper bound for the number of isolated individuals and another one for the final epidemic size under some pre–specified conditions.

Proposition 3 Let K_{Q} = \frac{θN}{μ} (1−e^{−μρ})t and K_{I} = \frac{βμN}{(μ − θ(1 − e^{−μρ}))(μ + γ + θ)}.

Then K_{Q} is an upper bound for the number of isolated individuals in the time interval [0, t] and K_{I} is an upper bound for the final epidemic size provided that

ρ < −\frac{1}{μ} ln (1 − \frac{μ}{θ}) .

4 Stability analysis of E^{∗} and the existence of Hopf bifurcation

In this section, we analyze the local stability of the endemic steady state for system (2) and look for the possibilities of stability switches. When analyzing the local stability of system (2), there are two cases to analyze (either ρ = 0 or ρ > 0). The case ρ = 0 implies that θ = 0. This corresponds to the case where there is no isolation and therefore, system (2) turns out to be the classical SIR model in this case and the detailed analysis of the classical SIR model can be found in [3, 19, 27]. Hence, we here focus on the case where ρ > 0.

Linearization of system (2) about the endemic steady state E^{∗} produces a variational system of the form

\frac{dx}{dt} = F(x, x_{ρ}) = Lx + Mx_{ρ},  (17)

where x and x_{ρ} denote the vectors of instantaneous and delayed variational variables, respectively, while L and M stand for the Jacobians given by

L = \frac{∂F(x, x_{ρ})}{∂x}|_{x=E_{1}} and M = \frac{∂F(x, x_{ρ})}{∂x_{ρ}} |_{x=E_{1}}. Determining the nature of

the roots of the characteristic equation of (17)

|λI − L − Me^{−λρ}| = 0,  (18)

is not trivial but is carried out below.

In addition, we show that Hopf bifurcation occurs at ρ = ρ_{0} > 0 whenever

λ = iω_{0} are roots of (18), where ω_{0} > 0 and \frac{d Re(λ)}{dρ}|_{ρ=ρ_{0}} ≠ 0 (see [17, 8]).

We now compute the variational system associated with (2) about E^{∗}. To this end, we assume S(t) = A(t) + S^{∗}, I(t) = B(t) + I^{∗}, Q(t) = C(t) + Q^{∗}.

On substituting in (8) and ignoring nonlinear terms, we obtain

A′(t) = (−a − μ)A(t) − bB(t) − cC(t),
B′(t) = aA(t) + (b − k)B(t) + cC(t),  (19)
C′(t) = θB(t) + dB(t − ρ) − μC(t),

where

a = β \frac{I^{∗}}{N − Q^{∗}} ,  b = β \frac{S^{∗}}{N − Q^{∗}} ,  c = β \frac{S^{∗}I^{∗}}{(N − Q^{∗})^{2}} ,
d = −θe^{−μρ},  and  k = θ + γ + μ.

System (19) can be written in the matrix form

\frac{dx}{dt} = Lx + Mx_{ρ},  (20)

where

x = [\begin{matrix}A \\ B \\ C \end{matrix}] ,  L = [\begin{matrix}−a − μ & −b & −c \\ a & b − k & c \\ 0 & θ & −μ \end{matrix}]  and  M = [\begin{matrix}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & d & 0 \end{matrix}].

The characteristics equation corresponding to system (20) is given by

|λI − L − Me^{−λρ}| = 0,  (21)

where | · | denotes matrix determinant. This characteristic equation could be written in the form

P(λ) + Q(λ, ρ)e^{−λρ} = 0,  (22)

where

P(λ) = a_{3}λ^{3} + a_{2}λ^{2} + a_{1}λ + a_{0},  Q(λ, ρ) = −cdλ − cdμ,  (23)

and,

a_{3} = 1,  a_{1} = ak + aμ + μ^{2} + 2kμ − 2bμ − θc,
a_{2} = a + k + 2μ − b,  a_{0} = (k − b)μ^{2} + akμ − μθc.

We will first establish the stability of the endemic equilibrium E^{∗} in the absence of delay (i.e., ρ = 0) and then consider the more general case ρ ≠ 0.

4.1 Analysis in the absence of delay (ρ = 0)

If we assume ρ = 0, then the characteristic equation (22) reads

λ^{3} + μ(R_{0} + 1)λ^{2} + μ [μ + (k + μ)(R_{0} − 1)] λ + μ^{2}k(R_{0} − 1) = 0. (24)

Moreover, the Hurwitz determinants are Δ_{1} = μ(R_{0} + 1), Δ_{2} = μ^{2}R_{0} [μ(R_{0} + 1) + k(R_{0} − 1)] and Δ_{3} = μ^{2}k(R_{0}−1)Δ_{2} which are all positive for R_{0} > 1. Thus, based on applying Descarte’s rule of signs and Routh–Hurwitz criterion, the roots of equation (24) will have negative real part for R_{0} > 1. Therefore, we show the following proposition.

Proposition 4 In the absence of delay (i.e., ρ = 0), the endemic equilibrium E^{∗} is locally asymptotically stable whenever it exists.

4.2 Analysis in the presence of delay (ρ ≠ 0)

Since e^{−μρ} = 1 − μρ + O ((μρ)^{2}), then (on assuming that μρ ≪ 1) we get

e^{−μρ} ∼ 1 − μρ.

Hence,

\frac{S^{∗}}{N} = \frac{1 − Θ̂μρ}{R_{0} − Θ̂μρ},  \frac{I^{∗}}{N} = \frac{μ}{k} \frac{R_{0} − 1}{R_{0} − Θ̂μρ},  \frac{Q^{∗}}{N} = \frac{Θ̂μρ(R_{0} − 1)}{R_{0} − Θ̂μρ}, (25)

and

a = \frac{μ(R_{0} − 1)}{1 − Θ̂μρ},  b = k,  c = \frac{a}{R_{0}} .  (26)

Thus,

a_{2} = 2μ + \frac{μ(R_{0} − 1)}{1 − Θ̂μρ},
a_{1} = μ^{2} + (k + μ − \frac{θ}{R_{0}} ) \frac{μ(R_{0} − 1)}{1 − Θ̂μρ},  (27)
a_{0} = (k − \frac{θ}{R_{0}} ) \frac{μ^{2}(R_{0} − 1)}{1 − Θ̂μρ}.

Also,

d = −θe^{−μρ} ∼ θμρ − θ.

Similarly, Q(λ, ρ) becomes

Q(λ, ρ) = b_{1}λ + b_{0},

where

b_{1} = θc − cθμρ = μ(1 − μρ)\frac{θ(R_{0} − 1)}{R_{0}(1 − Θ̂μρ)},  and
b_{0} = μθc − cθμ^{2}ρ = μb_{1}.  (28)

It is clear that the function P is a cubic polynomial, while the function Q is a linear polynomial in λ. Their coefficients a_{i} and b_{i} are functions in the delay parameter ρ. Now, we follow the same approach shown in [1] to check the stability switching possibility, as the delay parameter ρ is allowed to increase.

To this end, the characteristic equation (21) is rewritten in the form

D(λ, ρ) = 0,  (29)

where

D(λ, ρ) = P_{n}(λ, ρ) + Q_{m}(λ, ρ)e^{−λρ},  (30)

and where

P_{n}(λ, ρ) = ∑_{k=0}^{n} p_{k}(ρ)λ^{k},  Q_{m}(λ, ρ) = ∑_{k=0}^{m} q_{k}(ρ)λ^{k}.  (31)

with n = 3, m = 1. It is noteworthy that the coefficients p_{k}(ρ) = a_{k} and q_{k}(ρ) = b_{k} are rational functions that do not have any singularities for μρ ≪ 1. Thus, in (31), we have n, m ∈ ℕ^{∗}, n > m, and p_{k}(.), q_{k}(.) : ℝ^{+} → ℝ are continuous and differentiable functions of ρ such that

P_{n}(0, ρ)+Q_{m}(0, ρ) = p_{0}(ρ)+q_{0}(ρ) = \frac{μ^{2}(R_{0} − 1)}{1 − Θ̂μρ}(k − \frac{μρθ}{R_{0}})≠0, ∀ρ∈ℝ^{+}.

(32)

Hence, λ = 0 is not a characteristic root of (29). On the other hand, the functions P and Q are polynomials in λ. Therefore, they are analytic functions in a right half plane Re (λ) > −δ, δ > 0. Moreover, they are differentiable in ρ, and we need to verify the following conditions.

1) If λ = iω, ω ∈ ℝ, then P_{n}(iω, ρ) + Q_{m}(iω, ρ) ≠ 0, ρ ∈ ℝ.

2) lim sup \{|\frac{Q_{m}(λ, ρ)}{P_{n}(λ, ρ)}| : |λ| → ∞, ℜ(λ) ≥ 0\} < 1, for any ρ.

3) F (ω, ρ) := |P_{n}(iω, ρ)|^{2}−|Q_{n}(iω, ρ)|^{2} for each ρ has at most a finite number of real zeros.

4) Each positive root ω(ρ) of F (ω, ρ) = 0 is continuous and differentiable in ρ whenever it exists.

If these conditions hold, the following statements are true [1, 2].

(a) Suppose that the equation F (ω, ρ) = 0 has no positive zeros. Then if D(λ, ρ) is unstable at ρ = 0, it remains unstable for all ρ ≥ 0.

(b) Suppose that the equation F (ω, ρ) = 0 has at least one positive zero and that each positive zero is simple. As ρ increases, stability switches may occur. Moreover, there exists a positive number ρ^{∗} such that D(λ, ρ) is unstable for all ρ > ρ^{∗}. As ρ varies from 0 to ρ^{∗}, at most a finite number of stability switches may occur.

Now, on putting λ = iω, we get

P_{3}(iω, ρ) = (−a_{2}ω^{2} + a_{0}) + i(−ω^{3} + a_{1}ω), Q_{1}(iω, ρ) = b_{0} + ib_{1}ω. (33)

The condition 1) is obviously true because

P_{3}(iω, ρ) + Q_{1}(iω, ρ) = (b_{0} + a_{0} − ω^{2}a_{2}) + iω(b_{1} + a_{1} − ω^{2}) ≠ 0. (34)

For more explanation, the expression in (34) equals zero if and only if the expression (a_{2} − μ) [a_{0} + b_{0} + μ^{2}a_{2}] vanishes, which is impossible to occur. Thus,

(34) is true. Also, from (23) we have that Q_{m} and P_{n} are polynomials of different degrees in λ, with m < n. Thus,

|\frac{Q_{1}(λ, ρ)}{P_{3}(λ, ρ)}| → 0 as |λ| → +∞.

This proves condition 2). Now, let F be defined as in condition 3). Thus,

F (ω, ρ) = |P_{3}(iω, ρ)|^{2} − |Q_{1}(iω, ρ)|^{2},
= ω^{6}+(a_{2}^{2} − 2a_{1})ω^{4}+(a_{1}^{2} − b_{1}^{2} − 2a_{0}a_{2})ω^{2}+a_{0}^{2}−b_{0}^{2}, (35)
= 0.

Let z = ω^{2} and denote

p = a_{2}^{2} − 2a_{1},  q = a_{1}^{2} − b_{1}^{2} − 2a_{0}a_{2},  r = a_{0}^{2} − b_{0}^{2}.  (36)

Then, equation (35) becomes

z^{3} + pz^{2} + qz + r = 0,  (37)

where

p = 2μ^{2} + \frac{μ(R_{0} − 1)}{1 − Θ̂μρ} \{μ (\frac{R_{0} − 1}{1 − Θ̂μρ} + 2) − 2k (1 − \frac{Θ̂}{R_{0}} )\} , (38)

q = [μ^{2} − \frac{μ(R_{0} − 1)}{1 − Θ̂μρ} (k (1 − \frac{Θ̂}{R_{0}} ) − μ)]^{2} +

\frac{2μ^{3}k(R_{0} − 1)^{2}}{(1 − Θ̂μρ)^{2}} [1 − \frac{Θ̂}{R_{0}} − \frac{Θ̂}{R_{0}} \frac{θ(1 − μρ)^{2}}{2μR_{0}} ] ,  (39)
r = a_{0}^{2} − b_{0}^{2},
= \frac{μ^{4}(R_{0} − 1)^{2}}{(1 − Θ̂μρ)^{2}} [k − \frac{θ}{R_{0}} + \frac{(1 − μρ)θ}{R_{0}} ][k − \frac{θ}{R_{0}} − \frac{(1 − μρ)θ}{R_{0}} ],

= \frac{k^{2}μ^{4}(R_{0} − 1)^{2}}{(1 − Θ̂μρ)^{2}} [1 − \frac{μρΘ̂}{R_{0}}] [1 − \frac{Θ̂}{R_{0}} − \frac{(1 − μρ)Θ̂}{R_{0}}] .  (40)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 49–71, Jan–Jun 2020

Denote

H(z) = z^{3} + pz^{2} + qz + r.  (41)

Since H is a cubic polynomial in z, it has at most three roots. Also, since lim_{z→+∞} H(z) = +∞, we conclude that if r < 0, then equation (37) has at least

one (but no more than three) positive root(s). Hence, F (ω, ρ) = 0 has a finite number of real zeros and, therefore, condition 3) holds. Moreover, by the implicit function theorem, condition 4) is satisfied. Thus, F (ω, ρ) = 0 has a solution ω(ρ) > 0 if and only if r < 0.

Clearly, the first two terms in the right hand side of (40) are positive, while the last term could be positive or negative. Thus, the sign of the last term of (40) determines that of r. I.e., r < 0 if and only if

0 ≤ ρ < \frac{1}{μ} [2 − \frac{R_{0}}{Θ̂} ] , with \{1 ≤ R_{0} < 2Θ̂ and \frac{θ}{γ + μ} > 1\} . (42)

The conditions (42) say that for high enough isolation rate and for values of R_{0} slightly above 1 (e.g. 1 < R_{0} < 2), there is a range for the time delay ρ during which there is a positive solution for (37). However, if one of the inequalities (42) is broken, then r > 0.

On the other hand, equation (37) could be rewritten as

(z − \frac{−p}{2})^{2} + (q − \frac{p^{2}}{4} ) = \frac{−r}{z} .  (43)

We look for positive solutions of (43). The left hand side has graph of a vertical parabola with positive leading coefficient, which means that it opens up. Its vertex is at (−p/2, q − p^{2}/4). We have two cases.

Case 1: r < 0.

In this case, the right hand side of (43) is a monotonically decreasing function which is concave up. Its curve lies totally in the first quadrant. Moreover, it tends to ∞ as z tends to zero, while it tends to zero as z approaches ∞. Thus, the parabola will always cut this curve, either in one point or in three points representing the solutions of (37). Both curves will meet in three points if and only if the following set of inequalities is held

r < 0, p < 0, q > 0, p^{2} − 3q > 0.  (44)

Otherwise, for r < 0, there is a unique solution of (37). The case of having a unique positive solution corresponds to either of the following two sets of inequalities

{p > 0, r < 0} or {p < 0, q < 0, r < 0} .  (45)

Case 2: r > 0.

In this case, the right hand side of (43) is a monotonically increasing function which is concave down. Its curve lies totally in the fourth quadrant. It tends to −∞ as z tends to zero, while it tends to zero as z approaches ∞. Thus, the parabola will either cut this curve in two points or in no point. Thus, for r > 0, equation (37) will either have two positive solutions or will not have any positive solution. It has two positive solutions if and only if either of the following two inequality sets is held

{r > 0, p < 0, q > 0, p^{2} − 3q > 0} or {r > 0, q < 0, p > 0} . (46)

If neither of the above cases holds, then equation (37) doesn’t have a positive solution.

In the following proposition we summarise the above-mentioned conditions required for the existence of positive solutions of equation (37).

Proposition 5 Equation (37) has,

1. a unique positive solution if and only if the conditions (45) hold;

2. two positive solutions if and only if either of the two sets of inequalities in (46) holds;

3. three positive solutions if and only if all inequalities in (44) are satisfied;

4. no positive solution if neither of the above holds.

Algebraic manipulations aiming at combining and simplifying the conditions in (45), (46) and (44) would produce the sufficient conditions on model parameters showing the space of parameters on which at least one positive solution of (37) does exist. In case none of the conditions (45), (46) and (44) is satisfied, then there is no positive solution for (37) and therefore a stability switch cannot occur.

Let us now assume that λ = iω (with ω > 0) is a root of equation (22). Substituting and separating the real and imaginary parts result in

b_{1}ω cos (ωρ) − b_{0} sin (ωρ) = ω^{3}a_{3} − ωa_{1},
b_{0} cos (ωρ) + b_{1}ω sin (ωρ) = ω^{2}a_{2} − a_{0}.

Hence,

sin (ωρ) = \frac{(a_{2}ω^{2} − a_{0})b_{1}ω + (a_{1}ω − ω^{3})b_{0}}{b_{0}^{2} + b_{1}^{2}ω^{2}} ,
cos (ωρ) = \frac{(a_{2}ω^{2} − a_{0})b_{0} + (ω^{3} − a_{1}ω)b_{1}ω}{b_{0}^{2} + b_{1}^{2}ω^{2}} .  (47)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 27(1): 49–71, Jan–Jun 2020

On using (47) in the relation

sin^{2} (ωρ) + cos^{2} (ωρ) = 1,

we get

ω^{6} + pω^{4} + qω^{2} + r = 0,  (48)

which means that ω^{2} is a positive solution of (37), or equivalently (35). I.e., ω(ρ) is a real solution of F (ω, ρ) = |P_{3}(iω, ρ)|^{2} − |Q_{1}(iω, ρ)|^{2} = 0.

Now, we assume that z^{∗} is a positive solution of (37), where it satisfies the inequalities of the items of proposition 5. Moreover, we assume that

I = {ρ : ρ ≥ 0, and z^{∗} does exist} ,

which means that I ⊆ ℝ^{+} is the set of all time delays ρ for which ω(ρ) is a positive root of (35). Assume further that a value ρ ∈ I has been chosen such that ω(ρ) satisfies (35). Then, following [1], we can define the angle α(ρ) ∈ [0, 2π] such that sin α(ρ) and cos α(ρ) are given by the right hand sides of (47). I.e.,

sin α(ρ) = \frac{(a_{2}ω^{2} − a_{0})b_{1}ω + (a_{1}ω − ω^{3})b_{0}}{b_{0}^{2} + b_{1}^{2}ω^{2}} ,
cos α(ρ) = \frac{(a_{2}ω^{2} − a_{0})b_{0} + (ω^{3} − a_{1}ω)b_{1}ω}{b_{0}^{2} + b_{1}^{2}ω^{2}} .  (49)

Hence, a stability switch occurs, through the roots λ = ±iω(ρ), where ω(ρ) > 0 and is a solution of (35), at the values ρ

ρ_{n} = \frac{α(ρ) + 2nπ}{ω(ρ)} ,  (50)

for n ∈ ℕ^{∗} := {1, 2, 3, …}. Thus, for each n ∈ ℕ^{∗}, the relation (50) defines the maps ρ_{n} : I → ℝ^{+}, and the stability switch may occur only for the ρ values at which

ρ_{n} = ρ, for some n ∈ ℕ^{∗}.  (51)

Hence, the solution of (35) and (49) defines the maps (50), and the occurrence of stability switches takes place at the zeros of the functions

T_{n}(ρ) := ρ − ρ_{n}(ρ), n ∈ ℕ^{∗}.  (52)

Thus, based on applying [1, Thm 2.2 and Remark 2.1] of Beretta and Kuang, the statements of the following proposition are concluded.

Proposition 6 According to the above analysis, we conclude the following.

1. If none of the set of inequalities in (45), (46) and (44) is satisfied, then there is no positive solution of (37) and therefore no stability switch may occur (i.e., if the endemic equilibrium E^{∗} of (8) is stable for ρ = 0, see proposition 4), then it remains stable for all ρ > 0.

2. If model parameters (except ρ) are chosen such that conditions (45) are satisfied, then equation (37) will have a unique positive real solution (say ω_{1}) and therefore one stability switch will occur at the real root (say ρ_{1}) of T_{n}(ρ) = 0 for some n ∈ ℕ^{∗}. Thus, equation (22) will have a simple conjugate pair of purely imaginary eigenvalues ±iω_{1}(ρ_{1}) and hence E^{∗} will be asymptotically stable for ρ ∈ [0, ρ_{1}) and unstable for ρ > ρ_{1}.

3. If model parameters (except ρ) are chosen such that conditions (46) are held, then equation (37) will have two positive real solutions, say ω_{21} and ω_{22} with ω_{21} < ω_{22}. Hence, stability switches may depend on all real roots of T_{n}^{21}(ρ) = 0 (that corresponds to the root ω_{21}) and T_{n}^{22}(ρ) = 0 (corresponding to the root ω_{22}). If we assume ρ_{21} and ρ_{24} are positive roots for T_{0}^{22}(ρ) = 0, while ρ_{22} and ρ_{23} are positive roots for T_{0}^{21}(ρ) = 0, and no real roots for T_{0}^{22}(ρ) = 0 and T_{0}^{21}(ρ) = 0 when n>0, then by [1, Remark 4.2] we have ρ_{21} < ρ_{22} < ρ_{23} < ρ_{24}. Moreover, the endemic equilibrium E^{∗} will be asymptotically stable for ρ ∈ [0, ρ_{21}) ∪ (ρ_{22}, ρ_{23}) ∪ (ρ_{24}, ∞), while it will be unstable for ρ ∈ (ρ_{21}, ρ_{22}) ∪ (ρ_{23}, ρ_{24}).

4. If model parameters (except ρ) are chosen such that conditions (44) are satisfied, then equation (37) will have three positive solutions, say 0 < ω_{31} < ω_{32} < ω_{33}, and therefore stability switches may depend on all real roots of their corresponding sequences T_{n}^{31}(ρ) = 0, T_{n}^{32}(ρ) = 0 and T_{n}^{33}(ρ) = 0. Thus, a number of stability switches may happen.

5 Discussion

Disease patterns for human populations have been observed for some diseases over the past two centuries. It has been documented that infectious agents like measles or highly variable pathogens such as those responsible for the “flu” are linked to recurrent epidemic outbreaks. Hence, it is naturally reasonable to see that efforts to determine the mechanisms behind the recurrent epidemic outbreaks have been carried out for over a century with the aid of mathematical tools.

These mechanisms are the re-stocking of a large enough pool of susceptible individuals, a process directly tied in to the study of the role of age-structure on disease dynamics. Additional mechanisms include seasonality, delays from various natural epidemiological starts (asymptomatic or latent stages), behavioral changes in disease dynamics and delays driven by interventions (quarantine/isolation) [5, 4, 20, 19, 28, 9, 25, 24]. In this paper we have systematically derived the system of integro-differential equations model given by (7) and showed that system (7) is equivalent to the system of delay differential equations given by (8) for non-exponential waiting times in the Q class. In addition to showing that solutions with positive initial conditions are positive and bounded, we established sufficient conditions for the (i) uniqueness of the solutions, (ii) global stability of the disease free equilibrium when R_{0} ≤ 1, (iii) providing thresholds for stability switches. Moreover, during the the local stability analysis of the endemic equilibrium, we found a third order characteristic equation of system delay equation model (8). The analysis of the characteristic equation of system (8) shows that periodic solutions with sustained oscillation exist and stability switches occur when we have a fixed period of isolation time (mandatory isolation period).

Acknowledgements

Authors thank the reviewers as well as the editors for their effort.

References

[1] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay dependent parameters, SIAM J. Math. Anal. 33 (2002), no. 5, 1144–1165.

[2] F.G. Boese, Stability with respect to the delay: On a paper of K.L. Cooke and P. van den Driessche, J. Math. Anal. and Applications 228 (1998), no. 2, 293–321.

[3] F. Brauer, C. Castillo-Chávez, Mathematical models in population biology and epidemiology, Texts in Applied Mathematics 40, Springer-Verlag, New York, 2001.

[4] C. Castillo-Chavez, H.W. Hethcote, V. Andreasen, S.A. Levin, W.M. Liu, Cross-immunity in the dynamics of homogeneous and heterogeneous populations, in: T.G. Hallam, L.G. Gross & S. A. Levin (Eds.) Mathematical Ecology, Proceedings of the Autumn Course Research Seminars, (Trieste, 1986), World Sci. Publishing, Teaneck, NJ (1988), 303–316.

[5] C. Castillo-Chavez, H.W. Hethcote, V. Andreasen, S.A. Levin, W.M. Liu, Epidemiological models with age structure, proportionate mixing, and cross–immunity, J. Math. Biol. 27 (1989), no. 3, 233–258.

[6] G. Chowell, P.W. Fenimore, M.A. Castillo-Garsow, C. Castillo-Chavez, SARS outbreaks in Ontario, Hong Kong and Singapore: the role of diagnosis and isolation as a control mechanism, J. Theor. Biol. 224 (2003), no. 1, 1–8.

[7] R. Curtiss III, The impact of vaccines and vaccinations: challenges and opportunities for modelers, Math. Biosc. Eng. 8 (2011), no. 1, 77–93.

[8] O. Diekmann, S.A. van Gils, S.M. Verduyn Lunel, H.-O. Walther, Delay equations: Functional-, Complex- and Nonlinear Analysis, Applied Mathematical Sciences 110, Springer–Verlag, New York, 1995.

[9] M. Erdem, M. Safan, C. Castillo-Chavez, Mathematical models of influenza with imperfect quarantine, Bull. Math. Biol. 79 (2017), no. 7, 1612–1636.

[10] Z. Feng, W. Huang, C. Castillo-Chavez, On the role of variable latent periods in mathematical models for tuberculosis, J. Dynamics and Differential Equations 13 (2001), no. 2, 425–452.

[11] Z. Feng, H.R. Thieme, Recurrent outbreaks of childhood diseases revisited: the impact of isolation, Math. Biosci. 128 (1995), no. 1–2, 93–130.

[12] Z. Feng, H.R. Thieme, Endemic models with arbitrarily distributed periods of infection II: fast disease dynamics and permanent recovery, SIAM Journal on Applied Mathematics 61 (2000), no. 3, 983–1012.

[13] Z. Feng, H.R. Thieme, Endemic models with arbitrarily distributed periods of infection I: fundamental properties of the Model, SIAM Journal on Applied Mathematics 61 (2000), no. 3, 803–833.

[14] Z. Feng, Y. Yang, D. Xu, P. Zhang, M.M. McCauley, J.W. Glasser, Timely identification of optimal control strategies for emerging infectious diseases, J. Theor. Biol. 259 (2009), no. 1, 165–171.

[15] L. Gao, H. Hethcote, Simulations of rubella vaccination strategies in China, Math Biosc 202 (2006), no. 2, 371–385.

[16] G.F. Gensini, M.H. Yacoub, A.A. Conti, The concept of quarantine in history: from plague to SARS, Journal of Infection 49 (2004), no. 4, 257–261.

[17] J.K. Hale, S.M. Verduyn-Lunel, Introduction to Functional Differential Equations, Applied Mathematical Sciences 99, Springer-Verlag, New York, 1993.

[18] M. A. Herrera-Valdez, M. Cruz-Aponte, C. Castillo-Chavez, Multiple outbreaks for the same pandemic: local transportation and social distancing explain the different “waves” of A-H1N1pdm cases observed in Mexico during 2009, Math. Biosc. Eng. 8 (2011), no. 1, 21–48.

[19] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000), no. 4, 599–653.

[20] H.W. Hethcote, M. Zhien, L. Shengbing, Effects of quarantine in six endemic models for infectious diseases, Math. Biosci. 180 (2002), no. 1–2, 141–160.

[21] N.W. McLachlan, Modern operational calculus with applications in technical mathematics, Revised edition, Dover Publications, Inc., New York, 1962.

[22] A. Mubayi, C. Kribs-Zaleta, M. Martcheva, C. Castillo-Chavez, A costbased comparison of quarantine strategies for new emerging diseases, Math. Biosc. Eng. 7 (2010), no. 3, 687–717.

[23] M. Nuño, A mathematical model for the dynamics of influenza at the population and host level, Ph.D. Thesis, Cornell University, 2005.

[24] M. Safan, Mathematical analysis of an SIR respiratory infection model with sex and gender disparity: special reference to influenza A, Mathematical Biosciences and Engineering 16 (2019), no. 4, 2613–2649.

[25] M. Safan, M. Kretzschmar, K.P. Hadeler, Vaccination based control of infections in SIRS models with reinfection: special reference to pertussis, Journal of Mathematical Biology 67 (2013), no. 5, 1083–1110.

[26] E. Shim, Prioritization of delayed vaccination for pandemic influenza, Math. Biosc. Eng. 8 (2011), no. 1, 95–112.

[27] H.R. Thieme, Mathematics in population biology, Princeton Series in Theoretical and Computational Biology, Princeton University Press, Princeton, New Jersey, 2003.

[28] H.R. Thieme, The transition through stages with arbitrary length distributions, and applications in epidemics, in: C. Castillo-Chavez, S. Blower, P. van den Driessche, D. Kirschner, A.A. Yakubu (Eds.) Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory. (Minneapolis, MN, 1999), The IMA Volumes in Mathematics and its Applications, 126, Springer, New York, 2002, pp. 45–84.

[29] Y. Yang, Z. Feng, D. Xu, Analysis of a model with multiple infectious stages and arbitrarily distributed stage durations, Math. Model. Nat. Phenom. 3 (2008), no. 7, 180–193.