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

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

dengue model with early-life stage of vectors and age-structure within host

modelo de dengue con etapa temprana en vectores y estructura de edad

en el huésped

Fabio Sanchez^{∗}  Juan G. Calvo^{†}

Received: 7/Aug/2019; Revised: 26/Aug/2019; Accepted: 18/Sep/2019

∗University of Costa Rica, CIMPA-School of Mathematics, San José, Costa Rica. E-Mail: fabio.sanchez@ucr.ac.cr

†Misma dirección que/Same address as: F. Sánchez. E-Mail: juan.calvo@ucr.ac.cr

158 f. sanchez — j. g. calvo

Abstract

We construct an epidemic model for the transmission of dengue fever with an early-life stage in the vector dynamics and age-structure within hosts. The early-life stage of the vector is modeled via a general function that supports multiple vector densities. The basic reproductive number and vector demographic threshold are computed to study the local and global stability of the infection-free state. A numerical framework is implemented and simulations are performed.

Keywords: dengue fever; epidemic models; mathematical modeling; vectorborne diseases; age-structured model.

Resumen

Se construye un modelo epidémico para la dinámica de transmisión de la fiebre del dengue con etapa temprana de vida en los vectores y estructura de edad en las clases de los huéspedes. La etapa de vida temprana del vector, se modela a través de una función general que admite múltiples densidades de vectores. El número reproductivo básico y vector de umbral demográfico, se calculan para estudiar la estabilidad local y global del estado libre de infección. Se implementa además un marco numérico y se realizan simulaciones.

Palabras clave: dengue; modelos epidémicos; modelizaje matemático; enfermedades vectoriales; modelos con estructura de edad.

Mathematics Subject Classification: 93A30, 92B99, 37N25

1 Introduction

Dengue fever has been a burden to public health officials in the tropics and subtropics since the 20th century [2, 7]. Dengue virus, of the genus Flavivirus of the family Flaviviridae, is an infectious disease transmitted by the mosquitoes Aedes aegypti and Aedes albopictus [6]. There are four serotypes of the dengue virus, called DEN-1, DEN-2, DEN-3, and DEN-4. After infection of one serotype, the infected person acquires lifelong immunity for that specific serotype and shortterm immunity to other serotypes [2].

There are two stages in the transmission cycle of dengue that have been reported. There is an enzootic transmission cycle between primates mostly in forests, with transmission between vectors through feeding from infected animals [6]. These infected mosquitoes rarely wander far from the forest, so the infection to human populations comes from humans or livestock who visit a forest with presence of the virus, encounter an infected vector, become infected, and then infect the mosquitoes in their population center, which then can spread the disease to the rest of the population. In the case of rural, small populations, since the population rapidly gets saturated with the infection and subsequently immunized, the epidemic usually is short-termed. The other stage is between vectors and humans, where vectors bite an infected human and can potentially become infected. In this work we will focus on the interaction between vectors and humans.

The model in [11] focuses on the early-life stage of the vector and explores the effect of multiple vector densities on dengue outbreaks. Previous work on dengue models mostly focus on the adult vector-host interactions; see, e.g., [5, 3, 4, 12, 8, 9, 1, 10].

The model we present here incorporates age-structure within the host population, as well as the early-life stage of the vector as in [11]. The inclusion of age-structure in the human/host population can help to determine prevention and control strategies based on specific population age groups and other social factors inherent to a subgroup of the population at risk.

This article is organized as follows. In Section 2, we outline the mathematical model. In Section 3, we compute the vector demographic threshold, the basic reproductive number, and determine the stability of the system. Section 4 includes the numerical scheme and numerical simulations, confirming the theoretical results. Finally, in Section 5 we present some relevant conclusions and final thoughts.

2 Mathematical model

We consider a compartmental model with age-structure within the host population, with susceptible, infectious and recovered hosts, denoted by S_{h}(t, a), I_{h}(t, a) and R_{h}(t, a), respectively. Vectors are described by three states: E(t) (egg/larvae at time t), S_{v}(t) (number of non-infected vectors) and I_{v}(t) (number of infected vectors).

Hosts and vectors are coupled via a transmission process, where susceptible hosts can become infected at rate β(a)\frac{I_{v}(t)}{N_{v}(t)}, where β(a) represents the age-dependent contact rate (vector-human) and N_{v}(t) = S_{v}(t) + I_{v}(t) is the total number of vectors in the system. The number of new hosts coming into the system, Λ, is assumed to be constant. Infected individuals can recover at rate γ(a), and all hosts exit the system at rate μ_{h}(a).

We will restrict ourselves to the case of proportional mixing:

p(t, a) = \frac{c(a)n(t, a)}{∫_{0}^{∞} c(a)n(t, a) da},

with c(a) the age-specific per-capita contact/activity rate. We then define the force of infection

B(t) = ∫_{0}^{∞} \frac{I_{h}(t, a)}{n(t, a)}p(t, a)da,

where n(t, a) = S_{h}(t, a) + I_{h}(t, a) + R_{h}(t, a) is the population density and ∫_{0}^{∞} n(t, a)da is the total population.

In the vector classes, we consider that eggs enter the system at rate f(N_{v}), where f(N_{v}) is assumed to be a Kolmogorov type function, f(N_{v}) = N_{v}g(N_{v}), with g : ℝ^{+} → ℝ^{+} a differentiable function such that g(0) > 0, g(∞) = 0. They also exit the system at rate μ_{e} and become mosquitoes at rate δ. Mosquitoes become infectious at rate β_{v}B(t), where β_{v} is the transmission rate from humans to mosquitoes and B(t) is the force of infection. Mosquitoes also exit the system at rate μ_{v}. In our analysis, we assume that parameters δ, μ_{e}, μ_{v}, β_{v} are constant and β, μ_{h}, γ are continuous functions on age.

The model we just described is given by the system of differential equations:

\frac{dE}{dt} = f(N_{v}) − (δ + μ_{e})E,  (1a)
\frac{dS_{v}}{dt} = δE − β_{v}S_{v}B(t) − μ_{v}S_{v},  (1b)
\frac{dI_{v}}{dt} = β_{v}S_{v}B(t) − μ_{v}I_{v},  (1c)
(\frac{∂}{∂t} + \frac{∂}{∂a}) S_{h}(t, a) = −β(a)S_{h}(t, a)\frac{I_{v}}{N_{v}} − μ_{h}(a)S_{h}(t, a),  (1d)
(\frac{∂}{∂t} + \frac{∂}{∂a}) I_{h}(t, a) = β(a)S_{h}(t, a)\frac{I_{v}}{N_{v}} − (μ_{h}(a) + γ(a)) I_{h}(t, a), (1e)
(\frac{∂}{∂t} + \frac{∂}{∂a}) R_{h}(t, a) = γ(a)I_{h}(t, a) − μ_{h}(a)R_{h}(t, a),  (1f)

along with initial conditions given by

\begin{aligned}E(0) & = E_{0}, & S_{v}(0) & = S_{v_{0}}, & I_{v}(0) & = I_{v_{0}}, \\ S_{h}(t, 0) & = Λ, & I_{v}(t, 0) & = 0, & R_{h}(t, 0) & = 0, \\ S_{h}(0, a) & = S_{h_{0}}(a), & I_{h}(0, a) & = I_{h_{0}}(a), & R_{h}(0, a) & = R_{h_{0}}(a). \end{aligned}

The total host population satisfies

(\frac{∂}{∂t} + \frac{∂}{∂a}) n(t, a) = −μ_{h}(a)n(t, a),

and we then can compute explicitly that

n(t, a) = \{ \begin{matrix}n_{0}(a − t)\frac{F(a)}{F(a − t)} & if a ≥ t, \\ ΛF(a) & if a < t, \end{matrix}  (2)

where

F(a) = e^{− ∫_{0}^{a} μ_{h}(s) ds},

is the proportion of individuals that survive at age a. Therefore, we define

n^{∗}(a) := lim_{t→∞} n(t, a) = ΛF(a),
p_{∞}(a) := lim_{t→∞} p(t, a) = \frac{c(a)F(a)}{∫_{0}^{∞} c(a)F(a) da}.

3 Model analysis

In this section we explore the conditions for multiple vector demographic steady states and determine their stability. We also compute the basic reproductive number and analyze local and global stability for the solutions of system (1).

3.1 Vector demographic number, R_{v}

Since n(t, a) is given explicitly in (2), we first rescale variables

s_{h}(t, a) = \frac{S_{h}(t, a)}{n(t, a)}, i_{h}(t, a) = \frac{I_{h}(t, a)}{n(t, a)}, r_{h}(t, a) = \frac{R_{h}(t, a)}{n(t, a)},

to obtain the equivalent system

\frac{dE}{dt} = f(N_{v}) − (δ + μ_{e})E,  (3a)
\frac{dS_{v}}{dt} = δE − β_{v}B(t)S_{v} − μ_{v}S_{v},  (3b)
\frac{dI_{v}}{dt} = β_{v}B(t)S_{v} − μ_{v}I_{v},  (3c)
N_{v} = S_{v} + I_{v},  (3d)
(\frac{∂}{∂t} + \frac{∂}{∂a}) s_{h}(t, a) = −β(a)s_{h}(t, a)\frac{I_{v}}{N_{v}} ,  (3e)
(\frac{∂}{∂t} + \frac{∂}{∂a}) i_{h}(t, a) = β(a)s_{h}(t, a)\frac{I_{v}}{N_{v}} − γ(a)i_{h}(t, a),  (3f)
(\frac{∂}{∂t} + \frac{∂}{∂a}) r_{h}(t, a) = γ(a)i_{h}(t, a),  (3g)
B(t) = ∫_{0}^{∞} p(t, a)i_{h}(t, a)da.  (3h)

For a given equilibrium state (E^{∗}, S_{v}^{∗}, I_{v}^{∗}, N_{v}^{∗}, s_{h}^{∗}(a), i_{h}^{∗}(a), r_{h}^{∗}(a), B^{∗}) of system (3), we study its local stability by using the perturbations

E(t) = E^{∗} + e^{ψt}Ẽ,
S_{v}(t) = S_{v}^{∗} + e^{ψt}S̃_{v},
I_{v}(t) = I_{v}^{∗} + e^{ψt}Ĩ_{v},
N_{v}(t) = N_{v}^{∗} + e^{ψt}Ñ_{v},
s_{h}(t, a) = s_{h}^{∗}(a) + e^{ψt}s̃_{h}(a),
i_{h}(t, a) = i_{h}^{∗}(a) + e^{ψt}ĩ_{h}(a),
r_{h}(t, a) = r_{h}^{∗}(a) + e^{ψt}r̃_{h}(a),
B(t) = B^{∗} + e^{ψt}B̃,

where

B^{∗} = ∫_{0}^{∞} p_{∞}(a)i_{h}^{∗}(a) da, B̃ = ∫_{0}^{∞} p_{∞}(a)ĩ_{h}(a) da.

Linearization leads to the eigenvalue problem

ψẼ = f′(N_{v}^{∗})Ñ_{v} − (μ_{e} + δ) Ẽ,  (4a)
ψÑ_{v} = δẼ − μ_{v}Ñ_{v},  (4b)
ψĨ_{v} = β_{v} (S̃_{v}B^{∗} + S_{v}^{∗}B̃) − μ_{v}Ĩ_{v},  (4c)
ψS̃_{v} = δẼ − β_{v} (B^{∗}S̃_{v} + S_{v}^{∗}B̃) − μ_{v}S̃_{v},  (4d)
\frac{d}{da}s̃_{h}(a) + ψs̃_{h}(a) = −β(a) (\frac{I_{v}^{∗}}{N_{v^{∗}}} s̃_{h}(a) + \frac{Ĩ_{v}}{N_{v}^{∗}} s_{h}^{∗}(a) − \frac{I_{v}^{∗}}{N_{v}^{∗}} \frac{Ñ_{v}}{N_{v}^{∗}} s_{h}^{∗}(a)) ,  (4e)

\frac{d}{da}ĩ_{h}(a) + ψĩ_{h}(a) = β(a) (\frac{I_{v}^{∗}}{N_{v}^{∗}} s̃_{h}(a) + \frac{Ĩ_{v}}{N_{v}^{∗}} s_{h}^{∗}(a) − \frac{I_{v}^{∗}}{N_{v}^{∗}} \frac{Ñ_{v}}{N_{v}^{∗}} s_{h}^{∗}(a)) − γ(a)ĩ_{h}(a), (4f)
\frac{d}{da}r̃_{h}(a) + ψr̃_{h}(a) = γ(a)ĩ_{h}(a).  (4g)

For Ẽ, Ñ_{v} ≠ 0, equations (4a) and (4b) imply that

\frac{δf′(N_{v}^{∗})}{(μ_{e} + δ + ψ)(ψ + μ_{v})} = 1.

Let

H_{v}(ψ) = \frac{δf′(N_{v}^{∗})}{(μ_{e} + δ + ψ)(ψ + μ_{v})}.

We then define the demographic vector number

R_{v} = H_{v}(0) = \frac{f′(N_{v}^{∗})}{ϕ},

where

ϕ = \frac{(δ + μ_{e})μ_{v}}{δ} ,

represents the proportion of eggs that survive to the adult stage. Recall that the rate that eggs enters the system is given by f(N_{v}) = N_{v}g(N_{v}). We can then establish the following result:

Lemma 3.1 Suppose that the set g^{−1}(ϕ) = {N_{v} ∈ (0, +∞) : g(N_{v}) = ϕ} is non-empty. For each N_{v} ∈ g^{−1}(ϕ), there exists a positive vector state

(E^{∗}, N_{v}^{∗}) = (\frac{μ_{v}}{δ} N_{v}, N_{v}) ,  (5)

which is locally asymptotically stable if R_{v} < 1 and unstable otherwise.

Proof. We have that (E, N_{v}) satisfies the system

\frac{dE}{dt} = f(N_{v}) − (δ + μ_{e})E,
\frac{dN_{v}}{dt} = δE − μ_{v}N_{v},  (6)

with appropriate initial conditions. Therefore, fixed points satisfy

f(N_{v}^{∗}) = (δ + μ_{e})E^{∗},
δE^{∗} = μ_{v}N_{v}^{∗}.

Multiplying both equations and using the fact that f(N_{v}^{∗}) = N_{v}^{∗}g(N_{v}^{∗}), we then deduce that g(N_{v}^{∗}) = ϕ (for E^{∗}N_{v}^{∗} ≠ 0). Thus, for each N_{v} ∈ g^{−1}(ϕ) there exists the positive state given in (5).

Moreover, for a fixed state (5), the associated Jacobian to System (6) is given by

[ \begin{matrix}−(δ + μ_{e}) & f′(N_{v}^{∗}) \\ δ & −μ_{v} \end{matrix} ] ,

which eigenvalues are given by

\frac{1}{2} (−(δ + μ_{e} + μ_{v}) ± \sqrt{(δ + μ_{e} + μ_{v})^{2} − 4μ_{v}(μ_{e} + δ) (1 − R_{v})}) .

If R_{v} < 1, we then conclude that both eigenvalues have negative real part and (5) is locally stable.

Remark 3.2 Since we assume that f(N_{v}) = N_{v}g(N_{v}), it is straightforward to verify that

R_{v} = 1 + g′(N_{v}^{∗})\frac{N_{v}^{∗}}{ϕ} .

Thus, equilibrium points given by (5) are locally stable if g′(N_{v}^{∗}) < 0, and unstable otherwise; see Figure 1.

Figure 1: An example for g(N_{v}) as a function of N_{v} for which multiple steady states exist. The dashed line corresponds to the value of ϕ. Each intersection of both curves corresponds to an endemic state N_{v}^{∗} ∈ g^{−1}(ϕ). Black filled dots correspond to stable points since g′(N_{v}^{∗}) < 0 and open circles correspond to unstable fixed points; see Lemma 3.1 and Remark 3.2. In this case, N_{v}^{∗} = 0 is unstable; see Lemma 3.3.

Lemma 3.3 The vector-free state (E^{∗}, S_{v}^{∗}, I_{v}^{∗}) = (0, 0, 0) is locally asymptotically stable if g(0) < ϕ, and unstable otherwise.

Proof. Suppose that g(0) < ϕ. For N_{v}^{∗} = 0, R_{v} simplifies to R_{v} = \frac{g(0)}{ϕ} < 1.

Since H_{v} : [0, ∞) is a decreasing function of ψ, the equation H_{v}(ψ) = 1 can have only solutions with negative real part, and (E^{∗}, N_{v}^{∗}) is locally stable. The result then follows since S_{v}^{∗}, I_{v}^{∗} are non negative and N_{v}^{∗} = S_{v}^{∗}+I_{v}^{∗}. If g(0) > ϕ, then H_{v}(ψ) = 1 has one positive solution and the result holds.

3.2 Basic reproductive number, R_{0}

Consider now the disease-free state

(E^{∗}, S_{v}^{∗}, I_{v}^{∗}, N_{v}^{∗}, s_{h}^{∗}(a), i_{h}^{∗}(a), r_{h}^{∗}(a), B^{∗}) = (E^{∗}, S_{v}^{∗}, 0, S_{v}^{∗}, 1, 0, 0, 0),

for system (3). From (4c) we get

\frac{Ĩ_{v}}{N_{v}^{∗}} = \frac{β_{v}B̃}{ψ + μ_{v}} .

Substituting in (4f) and solving, we obtain that

ĩ(a) = \frac{β_{v}B̃}{ψ + μ_{v}} ∫_{0}^{a} β(τ)e^{− ∫_{τ}^{a}(ψ+γ(h)) dh} dτ.

Multiplying by p_{∞}(a) and integrating with respect to a, we deduce that

B̃ = \frac{β_{v}B̃}{ψ + μ_{v}} ∫_{0}^{∞} ∫_{0}^{a} p_{∞}(a)β(τ)e^{− ∫_{τ}^{a}(ψ+γ(h)) dh} dτ da.

For B̃ ≠ 0, we obtain that

\frac{β_{v}}{ψ + μ_{v}} ∫_{0}^{∞} ∫_{0}^{a} p_{∞}(a)β(τ)e^{− ∫_{τ}^{a}(ψ+γ(h)) dh} dτ da = 1.

Let

G(ψ) = \frac{β_{v}}{ψ + μ_{v}} ∫_{0}^{∞} ∫_{0}^{a} p_{∞}(a)β(τ)e^{− ∫_{τ}^{a}(ψ+γ(h)) dh} dτ da.

We then define the basic reproductive number

R_{0} = G(0) = \frac{β_{v}}{μ_{v}} ∫_{0}^{∞} ∫_{0}^{a} p_{∞}(a)β(τ)e^{− ∫_{τ}^{a} γ(h) dh} dτ da.  (7)

In the particular case of constant parameters, it reduces to

R_{0} = \frac{β_{v}β}{μ_{v}(μ_{h} + γ)}.

We then have the following results: Theorem 3.4 Assume that R_{0} < 1. Then, the disease-free solution of system (3) is globally asymptotically stable. Proof. From (3f) we have that

i_{h}(t, a) = ∫_{0}^{a} β(τ)e^{− ∫_{τ}^{a} γ(h) dh}s_{h} (τ + t − a, τ) \frac{I_{v}(τ + t − a)}{N_{v}(τ + t − a} dτ

for t > a. Multiplying by p(t, a) and integrating with respect to a we get

B(t) = ∫_{0}^{∞} ∫_{0}^{a} β(τ)p(t, a)e^{− ∫_{τ}^{a} γ(h) dh}s_{h}(τ + t − a, τ)\frac{I_{v}(τ + t − a)}{N_{v}(τ + t − a)} dτ.

Since s(t, a) ≤ 1,

B(t) ≤ ∫_{0}^{∞} ∫_{0}^{a} β(τ)p(t, a)e^{− ∫_{τ}^{a} γ(h) dh} \frac{I_{v}(τ + t − a)}{N_{v}(τ + t − a)} dτ,

and therefore

B^{∗} ≤ \frac{I_{v}^{∗}}{N_{v}^{∗}} ∫_{0}^{∞} ∫_{0}^{a} β(τ)p_{∞}(a)e^{− ∫_{τ}^{a} γ(h) dh} dτ.  (8)

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

From (3c), if S_{v}^{∗} = 0 then I_{v}^{∗} = 0. Otherwise, it holds that

B^{∗} = \frac{μ_{v}}{β_{v}} \frac{I_{v}^{∗}}{S_{v}^{∗}} .  (9)

Combining (8), (9), and using (7), we obtain that

0 ≤ \frac{I_{v}^{∗}}{N_{v}^{∗}} ≤ \frac{I_{v}^{∗}}{S_{v}^{∗}} ≤ \frac{I_{v}^{∗}}{N_{v}^{∗}} R_{0},

since S_{v}^{∗} ≤ N_{v}^{∗}. By assumption, R_{0} < 1, and therefore I_{v}^{∗} = 0. Hence, B^{∗} = 0 and i_{h}^{∗}(a) = 0.

Theorem 3.5 Assume that there exists N ∈ {N ∈ (0, +∞) : g(N) = ϕ} with g′(N) < 0. If R_{0} > 1, there exists one endemic non-uniform stable steady state for system (3).

Proof. It is straightforward to verify that there exists the endemic state

N_{v}^{∗} = N,
E^{∗} = N \frac{μ_{v}}{δ} ,
S_{v}^{∗} = N \frac{μ_{v}}{(μ_{v} + β_{v}B^{∗})},
I_{v}^{∗} = N \frac{β_{v}B^{∗}}{(μ_{v} + β_{v}B^{∗})}.

By hypothesis and Remark 3.2, it holds that (E^{∗}, N_{v}^{∗}) is a positive local stable fixed point. We will prove that B^{∗} > 0 which implies that I_{v}^{∗} > 0.

A non-uniform steady state for system (3) is a solution of the nonlinear system

\frac{I_{v}^{∗}}{N_{v}^{∗}} = \frac{β_{v}B^{∗}}{μ_{v} + β_{v}B^{∗}} ,
\frac{ds_{h}^{∗}(a)}{da} = −β(a)s_{h}^{∗}(a)\frac{I_{v}^{∗}}{N_{v}^{∗}} ,
\frac{di_{h}^{∗}(a)}{da} = β(a)s_{h}^{∗}(a)\frac{I_{v}^{∗}}{N_{v}^{∗}} − γ(a)i_{h}^{∗}(a),  (10)
\frac{dr_{h}^{∗}(a)}{da} = γ(a)i_{h}^{∗}(a),
B^{∗} = ∫_{0}^{∞} p_{∞}(a)i_{h}^{∗}(a) da,

for a > 0, with initial conditions

s_{h}^{∗}(0) = 1, i_{h}^{∗}(0) = 0, r_{h}^{∗}(0) = 0.

Consider the linear system of equations with parameter B given by

\frac{ds_{h_{B}(a)}^{∗}}{da} = −β(a)s_{h}^{∗}(a)\frac{β_{v}B}{(μ_{v} + β_{v}B)},
\frac{di_{h_{B}(a)}^{∗}}{da} = β(a)s_{h}^{∗}(a)\frac{β_{v}B}{(μ_{v} + β_{v}B)} − γ(a)i_{h}^{∗}(a),  (11)
\frac{dr_{h_{B}(a)}^{∗}}{da} = γ(a)i_{h}^{∗}(a),

for a > 0, with initial conditions

s_{h_{B}}^{∗}(0) = 1, i_{h_{B}}^{∗}(0) = 0, r_{h_{B}}^{∗}(0) = 0.

Given the solution (s_{B}^{∗}(a), i_{B}^{∗}(a), r_{B}^{∗}(a)) of system (11), define

H(B) = ∫_{0}^{∞} i_{B}^{∗}(a)p_{∞}(a) da.

It holds that (s_{B}^{∗}(a), i_{B}^{∗}(a), r_{B}^{∗}(a)) satisfies system (10) if and only if B is a fixed point of H; i.e., H(B) = B. Moreover, if B = 0 then i_{B}^{∗}(a) = 0 and H(0) = 0. Thus, in order to guarantee existence of at least one non-trivial solution to (10), it is just necessary to prove that H(B) has a positive fixed point.

Solving for s_{h_{B}}^{∗}(a) and i_{h_{B}}^{∗}(a), it follows that

s_{h_{B}}^{∗}(a) = e^{\frac{−βvB}{μv+βvB} ∫_{0}^{a} β(h) dh},
i_{h_{B}}^{∗}(a) = \frac{β_{v}B}{μ_{v} + β_{v}B} ∫_{0}^{a} e^{− ∫_{τ}^{a} γ(h) dh}β(τ)s_{h_{B}}^{∗}(τ) dτ.

Define

G(B) := \frac{H(B)}{B}, for B ≠ 0.

The function G : [0, 1] → ℝ is continuous by defining G(0) = H′(0). We have that

G(0) = \frac{β_{v}}{μ_{v}} ∫_{0}^{∞} ∫_{0}^{a} e^{− ∫_{τ}^{a} γ(h) dh}β(τ)p_{∞}(a) dτ da = R_{0} > 1,

and G(1) = H(1) < 1. Therefore, there exists B^{∗} ∈ (0, 1) such that G(B^{∗})=1, i.e., H(B^{∗}) = B^{∗} and there exists at least one endemic state for i_{h}^{∗}(a). Moreover, B^{∗} is unique since G(B) is strictly decreasing. In particular, B^{∗} > 0 implies that I_{v}^{∗} > 0, reaching an endemic steady state, both for vectors and humans.

4 Numerical implementation

For simplicity, we discretize the system of partial differential equations (3) with a first-order upwind finite difference scheme. We approximate the solution on the physical domain of interest given by the rectangle {(t, a) : t ∈ [0, T], a ∈ [0, A]}. We divide the intervals [0, T] and [0, A] in N_{T} and N_{A} subintervals, respectively, and consider the grid given by the nodes

(t_{j}, a_{k}) = (jΔt, kΔa) ,

for j ∈ {0, 1, …, N_{T} }, k ∈ {0, 1, …, N_{A}}, where

Δt := \frac{T}{N_{T}} , Δa := \frac{A}{N_{A}} ,

are the corresponding step sizes. For any function x and a grid point (t_{j}, a_{k}), we denote the approximation of x(t_{j}, a_{k}) by x_{k}^{j}. If the function depends only on age or time, it is denoted simply by x_{k} or x^{j}. We approximate the force of infection B(t_{j}) by the composite trapezium rule

B^{j} = Δa ( ∑_{k=1}^{N_{A}−1} p_{j}^{k} (i_{h})_{j}^{k} da + \frac{1}{2}p_{j}^{N_{A}}(i_{h})_{j}^{N_{A}}) .

We then have the explicit scheme given by the equations

(N_{v})^{j} = (S_{v})^{j} + (I_{v})^{j},
E^{j+1} = E^{j} + Δt [f ((N_{v})^{j}) − (δ + μ_{e})E^{j}] ,
(S_{v})^{j+1} = E^{j} + Δt [δE^{j} − β_{v}B^{j}(S_{v})^{j} − μ_{v}(S_{v})^{j}] ,
(I_{v})^{j+1} = (I_{v})^{j} + Δt [β_{v}B^{j}(S_{v})^{j} − μ_{v}(I_{v})^{j}] ,
(s_{h})_{k}^{j+1} = (s_{h})_{k}^{j} + Δt [−β_{k}(s_{h})_{k}^{j} \frac{(I_{v})^{j}}{(N_{v})^{j}} − \frac{(s_{h})_{k}^{j} − (s_{h})_{k−1}^{j}}{Δa} ] ,

(i_{h})_{k}^{j+1} = (i_{h})_{k}^{j} + Δt [β_{k}(s_{h})_{k}^{j} \frac{(I_{v})^{j}}{(N_{v})^{j}} − γ_{k}(i_{h})_{k}^{j} − \frac{(i_{h})_{k}^{j} − (i_{h})_{k−1}^{j}}{Δa} ] ,

(r_{h})_{k}^{j+1} = (r_{h})_{k}^{j} + Δt [γ_{k}(i_{h})_{k}^{j} − \frac{(r_{h})_{k}^{j} − (r_{h})_{k−1}^{j}}{Δa} ] ,

for 1 ≤ k ≤ N_{A} and 0 ≤ j ≤ N_{T} − 1. Thus, given the initial conditions E^{0}, (S_{v})^{0}, (I_{v})^{0}, s_{0}^{j}, s_{k}^{0}, i_{0}^{j}, i_{k}^{0}, r_{0}^{j}, r_{k}^{0}, we can compute the values of the unknowns on the grid points in successive time steps.

We present different scenarios where we confirm the results from Theorems 3.4 and 3.5. For the age dependent parameters, we show the distributions used in Figure 2. We consider three different functions f(N_{v}) = N_{v}g(N_{v}) in the following sections: a logistic-type function, a case with multiple vector steadystates, and a seasonal example.

Figure 2: (Left) transmission rate β(a), (middle) recovery rate γ(a) and (right) mortality rate μ(a), as functions of age.

4.1 Logistic growth

We first consider

g(N_{v}) = r (1 − \frac{N_{v}}{N_{max}} ) ,

for given constants r (mosquito growth rate) and N_{max} (maximum number of mosquitoes that the system can hold); see Figure 3. Recall that for a positive steady-state on vectors we require ϕ = g(N_{v}). In this case:

1. If ϕ > r, there exists only the trivial state N_{v}^{∗} = 0, which is stable since

R_{v} = \frac{g(0)}{ϕ} < 1; see Figure 3a.

2. For ϕ < r, besides the unstable state N_{v}^{∗} = 0, we have the non-zero state

N_{v}^{∗} = N_{max} (1 − \frac{ϕ}{r}) ,

which is stable since g′(N_{v}^{∗}) < 0; see Remark 3.2 and Figure 3b.

(a) N_{v}^{∗} = 0 is the only (stable) fixed point.

(b) Existence of a positive local stable point N_{v}^{∗}.

Figure 3: g(N_{v}) (solid line) as a function of the number of vectors N_{v}, considered in Section 4.1. The dashed line corresponds to the value ϕ. Solutions to the equation g(N) = ϕ correspond to non-trivial steady states for N_{v}.

Example 4.1 We first confirm that the condition R_{0} > 1 could lead to an endemic state, as long as there is a stable positive steady-state for vectors. We consider a set of parameters for which ϕ ≈ 3.39 and R_{0} ≈ 1.60. First, if r = 0.20, the only stable fixed point is (E^{∗}, S_{v}^{∗}, I_{v}^{∗}) = (0, 0, 0), for which i^{∗}(a) = 0; see Figures 4a, 4b. Second, if r = 5, we have the stable fixed point (E^{∗}, S_{v}^{∗}, I_{v}^{∗}) ≈ (9633, 7505, 88824). In this case, we have an endemic state as shown in Figures 4c, 4d, according to Theorem 3.5.

(a) i_{h}(t, a).  (b) Vectors.

(c) i_{h}(t, a).  (d) Vectors.

Figure 4: Solutions for the (left) infected class and (right) vectors when R_{0} > 1 for (top) ϕ > r = 0.20 and (bottom) ϕ < r = 5; see Example 4.1.

Example 4.2 In this example we confirm that the condition R_{0} < 1 is sufficient to guarantee a disease-free steady state. We take r = 5 and reduce β such that R_{0} < 1. The infected class reaches a disease-free state as shown in Figure 5a, according to Theorem 3.4. Even though there exists a positive state for vectors (E^{∗}, N_{v}^{∗}) as shown in Figure 5b, we observe that I_{v}^{∗} = 0.

(a) i(t, a).  (b) Vectors.

Figure 5: When R_{0} < 1, the disease-free state is stable, even though there is a positive steady state for vectors (R_{v} > 1); see Example 4.2.

Example 4.3 We now consider the case R_{0} > 1 with initial conditions E_{0} = 0, S_{v_{0}} = 10, I_{v_{0}} = 1, i_{h}(0, a) = r_{h}(0, a) = 0 (no infected or immune humans at time t = 0); see results for the infected class in Figure 6. It is clear that R_{0} > 1 guarantees an endemic state as long as vectors can survive.

(a) i_{h}(t, a).  (b) Vectors.

Figure 6: When R_{0} > 1, if there exists a positive state for vectors (R_{v} > 1) we can observe an endemic state on humans; see Example 4.3.

4.2 Multiple vector demographic states

In a second set of experiments we use

g(N_{v}) = re^{−N_{v}/c_{1}}(sin(c_{2}N_{v}) + 1),  (12)

for given constants r, c_{1}, c_{2}. Here, r is the vector per-capita fertility rate, c_{1} is a form of vector control and c_{2} represents the variations in vector densities; for a particular choice of parameters see Figure 1. Equation (12) represents the different growth rates of vectors for the wet and dry seasons. In this way, we simulate variations based on vector control efforts, obtaining multiple vector demographic steady states for (E^{∗}, N_{v}^{∗}).

For the particular choice of parameters we have used, we obtain eight positive fixed points. Numerically we confirm that four of them are locally stable.

Example 4.4 We first confirm the result proved in Lemma 3.3. We observe that if r < ϕ, the infection-free steady state is stable and unstable otherwise; results for (E_{0}, S_{v0}, I_{v0}) = (10, 10, 10) are shown in Figure 7. We then obtain different solutions for different initial conditions for the vector classes for which different positive steady-states are reached; see Figure 8. Despite having multiple vector densities the outbreaks are similar in severity. This implies that even when vector density is low an outbreak is possible.

(a) Vector growth rate (r = 3) is less than ϕ.

(b) Vector growth rate (r = 4) is greater than ϕ.

Figure 7: Different steady state distributions of infected individuals i_{h}^{∗}(t, a) with initial conditions (E_{0}, N_{v0}) = (10, 20). If r < ϕ, the infection-free steady state is stable and unstable otherwise; see Example 4.4.

(a) (E_{0}, N_{v0}) = (1525, 27460).  (b) Vectors.

(c) (E_{0}, N_{v0}) = (200000, 100000).  (d) Vectors.

Figure 8: (Left) infected individuals i_{h}(t, a) and (right) vector solutions with different initial conditions. For this choice of parameters, R_{0} > 1, R_{v} > 1 and several vector steady states exist; see Example 4.4.

Example 4.5 Similarly as Example 4.2, we confirm that R_{0} < 1 is sufficient for the disease to die out, even in the presence of a positive population of vectors. In this case, (E^{∗}, N_{v}^{∗}) = (7430, 74305), but I_{v}^{∗} = 0; see Figure 9. Even when the demographic vector number is bigger than one the disease can be under control when R_{0} < 1.

Figure 9: R_{0} < 1 is sufficient to guarantee that (left) i_{h}(a) = 0 and (right) I_{v}^{∗} = 0, even though S_{v}^{∗} > 0, E^{∗} > 0; see Example 4.5.

Example 4.6 Similarly as Example 4.3, we confirm that R_{0} > 1 implies the existence of an endemic state, as long as a positive equilibrium state exists for the vectors to survive; see Figure 10.

Figure 10: (Left) i_{h}(t, a) and (right) (E(t), S_{v}(t), I_{v}(t)) for R_{0} > 1; see Example 4.6. Initially all humans are susceptible and (E_{0}, S_{v0}, I_{v0}) = (0, 10, 1).

4.3 Effect of seasonality on dengue dynamics

In most places where dengue is endemic, seasonal variations in vector populations play a major role in disease transmission. Moreover, it determines the distribution of resources allocated for preventive/control measures. Typically, dengue incidence is correlated with the rainy season. The importance of understanding seasonal variations per location could potentially help public health officials to allocate resources, as well as having better preventive/control measures to reduce dengue incidence (focused primarily towards the reduction of vector breeding sites).

Figure 11: Solution for i(t, a) when seasonal parameters are considered, for (left) β_{0} = 0.09 and (right) β_{0} = 0.19; see Section 4.3.

We include some numerical results where g(N_{v}), β_{v} and δ depend periodically on time, simulating high and low seasons in the dynamics of vectors. Here we use parameters for the vector classes as in [11]. We consider a population with only susceptible humans. In the vector classes, we include one infected vector in order to observe the propagation of the disease; see results in Figure 11. For different values of the transmission rate (β(a)) the infected host distribution distinctly affects the younger and senior age groups.

5 Discussion

We have constructed a model with age-structure within host and early-life stage of the vector with a general function f(N_{v}) that represents the new vectors from the egg/larvae stage in the vector system. This gives the possibility of multiple demographic steady states for the vector population and its stability depends on the vector demographic number, R_{v}. The local stability of the vector-free state when R_{v} < 1 was established. The basic reproductive number was computed and the local and global asymptotic stability of the disease-free equilibrium was determined when R_{0} < 1.

When R_{0} >1 and we have a stable vector demographic steady state (R_{v} >1), the disease is then endemic. Control measures on the early-life stage of the vector can guarantee an adult vector-free state and hence, the disease dies out. Vector control measures such that ϕ > max g(N_{v}) implies that vectors will die out independently on the value of R_{0}.

There are important public health implications when we are able to include host age distribution, which can determine better strategies for hospitalized individuals. Furthermore, control measures on the early-life stage of the vector can effectively change the landscape on how public health officials lead prevention efforts before the onset of a dengue outbreak.

Acknowledgements

We thank the Research Center in Pure and Applied Mathematics and the Mathematics Department at Universidad de Costa Rica for their support during the preparation of this manuscript. The authors gratefully acknowledge institutional support for project B8747 from an UCREA grant from the Vice Rectory for Research at Universidad de Costa Rica.

References

[1] F. Brauer, C. Castillo-Chavez, A. Mubayi, S. Towers, Some models for epidemics of vector-transmitted diseases, Infect. Dis. Model. 1(2016), no.

1, 79–87. doi: 10.1016/ j.idm. 2016. 08.001

[2] Centers for Disease Control and Prevention, Dengue, https:// www.

cdc. gov/ Dengue/

[3] L. Esteva, C. Vargas, Analysis of a dengue disease transmission model, Math. Biosci. 150 (1998), no. 2, 131–151. doi: 10.1016/

S0025-5564(98)10003-2

[4] L. Esteva, C. Vargas, A model for dengue disease with variable human population, J. Math. Biol. 38(1999), no. 3, 220–240. doi: 10.1007/

s002850050147

[5] Z. Feng, J.X. Velasco-Hernández, Competitive exclusion in a vector-host model for the dengue fever, J. Math. Biol. 35(1997), no. 5, 523–544. doi:

10.1007/ s002850050064

[6] D.J. Gubler, Resurgent vector-borne diseases as a global health problem, Emerging Infect. Dis. 4(1998), no. 3, 442–450. doi: 10.3201/

eid0403. 980326

[7] E. Harris, E. Videa, L. Pérez, E. Sandoval, Y. Téllez, M.L. Pérez, R. Cuadra, J. Rocha, W. Idiaquez, R.E. Alonso, M.A. Delgado, L.A. Campo, F. Acevedo, A. Gonzalez, J.J. Amador, A. Balmaseda, Clinical, epidemiologic, and virologic features of dengue in the 1998 epidemic in Nicaragua, Am. J. Trop. Med. Hyg. 63(2000), no. 1-2, 5–11. doi: 10.4269/ ajtmh.

2000. 63.5

[8] C.A. Manore, K.S. Hickmann, S. Xu, H.J. Wearing, J.M. Hyman, Comparing dengue and chikungunya emergence and endemic transmission in A. aegypti and A. albopictus, J. Theor. Biol. 356(2014), 174–191. doi: 10.

1016/ j.jtbi. 2014. 04.033

[9] D. Murillo, S. Holechek, A. Murillo, F. Sanchez, C. Castillo-Chavez, Vertical transmission in a two-strain model of dengue fever, Letters in Biomathematics 1(2014), no. 2, 249–271. doi: 10.1080/ 23737867. 2014.

11414484

[10] F. Sanchez, L. Barboza, D. Burton, A. Cintrón-Arias, Comparative analysis of dengue versus chikungunya outbreaks in Costa Rica, Journal Ricerche di Matematica 67(2018), no. 1, 163–174. doi: 10.1007/

s11587-018-0362-3

[11] F. Sanchez, M. Engman, L. Harrington, C. Castillo-Chavez, Models for dengue transmission and control, in: A. Gumel, C. Castillo-Chavez, D.P. Clemence, & R.E. Mickens (Eds.) Mathematical studies on human disease dynamics. Emerging paradigms and challenges (Snowbird, Utah, 2005), Contemp. Math. 410(2006), 311–326. doi: 10.1090/ conm/

410/ 07734

[12] F. Sanchez, D. Murillo, C. Castillo-Chavez, Change in host behavior and its impact on the transmission dynamics of dengue, in: R.P. Mondani (Ed.) BIOMAT 2011, World Scientific, Singapore, 2012, pp. 191–203. doi: 10.

1142/ 9789814397711_0013