Dynamics of a dengue disease transmission model with two-stage structure in the human population

1.   Introduction

In recent years, dengue has become a problem for health systems at the international level. According to the World Health Organization (WHO), approximately between 50–100 million cases of dengue occur annually. Around 3.9 million people in approximately 128 countries are at risk of suffering from this disease because they live in regions where it is endemic. These countries are mainly located in regions of Africa, Asia, Oceania, and America where the geographic, demographic, and epidemiological conditions are conducive to the interaction between the vector, the virus, and the host. Dengue is transmitted by the bite of female mosquitoes, mainly of the species Aedes aegypti, after they become infected and after the minimum extrinsic incubation period has elapsed, which is 8 to 12 days before becoming infectious. The virus is injected into the human bloodstream by the mosquito bite and begins an incubation period for the disease that lasts from 4 to 10 days. After this time has elapsed, the infectious period of the host begins, starting the infectious process again ^[1].

For the study of viral dynamics, the mechanisms involved in the transmission of the pathogen between humans and vectors and the factors that influence the incidence of this disease, different mathematical models have been proposed in the language of ordinary differential equations (ODE). The results of such studies can aid the interventions by the health authorities.

In 1998, Esteva and Vargas ^[2] developed the first mathematical SIR-SI model for the dynamics of transmission of dengue fever that considers a constant human population.

Subsequently, Esteva and Vargas ^[3] formulated a system of differential equations that models the dynamics of transmission of dengue, where they consider a vertical and mechanical transmission in the vector population. Garba et al. ^[4] developed and rigorously analyzed a deterministic model for the dynamics of the transmission of dengue with or without imperfect vaccination, which exhibits a backward bifurcation.

Ghosh et al. ^[5] analyzed a compartmental model where the vector population is divided into three compartments: susceptible, exposed, and infected.

The human population is divided into seven compartments, within which two correspond to susceptible humans: those with low and high risk of getting infected. They assumed that the individuals in the low risk susceptible class have a lower chance of acquiring the infection as compared to the people in the high risk susceptible class.

For more details on the models in dengue, consult the systematic review conducted by Andraud et al. ^[6] in 2012, whose aim was to review the deterministic models of dengue transmission, and who identified 42 published models of interest, categorized into single serotype and multiserotype models.

Abidemi et al. ^[7] studied a two-strain compartmental dengue model that incorporates two control measures: Dengvaxia vaccine and insecticide. The same authors have studied other dengue model with multiple controls ^[8]. Asamoah et al. ^[9] analyzed a dengue transmission model in the presence of partial immunity and asymptomatic individuals with time-dependent controls.

They imposed time-dependent controls on their model: treated bednet, vaccination, treatment (prophylactics), and insecticides. Abidemi et al. ^[10] presented and analyzed a compartmental mathematical model that describes the dynamics of the transmission and control of dengue incorporating the effects of constant control rates. Abidemi et al. ^[11] developed and analyzed a model that captures the subpopulation of symptomatic infected human with severe dengue symptoms and Wolbachia-infected mosquito population.

There have been many studies of the incidence of dengue in children and adults. According to these studies, the most susceptible age groups vary, depending on the region. Siqueira Jr. et al. ^[12] studied dengue hemorrhagic fever (DHF) in Brazil and found that it primarily affected adults but noted an increase in cases in younger age groups. Huy et al. ^[13] performed an analysis of dengue data reported in Cambodia between the years 2002–2008, where they found that incidences by age were higher in infants under 1 year of age and children from 4 to 6 years old. In Thailand, dengue is more common in older children ^[14,15], while in Viet Nam the incidence is highest in children aged 6 to 10 years ^[16,17]. In Malaysia and Singapore it is more common in 18 year olds ^[18,19].

In 2016, Alera et al. ^[20] carried out a prospective longitudinal cohort study to determine the incidence of infection by the dengue virus in adults and children in the Philippines, where dengue is hyperendemic. They demonstrated that symptomatic dengue is primarily a pediatric disease in hyperendemic areas with a high force of infection.

In the reviewed literature on dengue models using ordinary differential equations, we identified some models that consider children and adults as different subpopulations. Aldila et al. ^[21] analyze an optimal control problem for a two-age-classes model, incorporating the implementation of mosquito repellent. Supriatna et al. ^[22] proposed a two-age-classes model with vaccination, because vaccination is usually concentrated in children. Chamnan et al. ^[23] developed a three-age-classes model with vaccination. The total population was subdivided into three groups: subjects less than 10 years old, between 10-44 years old, and those greater than 44 years old. This model ignores the terms of the transition from one age class to the next age class; this technical detail makes the constructions of the Lyapunov functions easier.

Although these dengue models take age groups into consideration, they have neglected the fact that adults are less susceptible to mosquito bites because they generally spend less time outdoors than children ^[24] and children are dependent upon adults for taking preventative measures. Translated into a measure of risk, we have that the relative risk of mosquito bites in adults versus children is less than unity. For this reason, we propose to analyze a two-stage-structure model that takes into consideration the differences in the risk of mosquito bites in children and adults.

This paper is organized as follows: in Section 2, we will propose a host–vector model with two-stage structure for dengue. In Section 3, we will compute the equilibria and basic reproductive number. In Section 4, we will study the local and global stability of disease-free equilibrium. In Section 5, we will determine the conditions for global asymptotic stability of the endemic equilibrium. In Section 6, we will present the results of numerical simulations for each of the scenarios of an outbreak of dengue. In Section 7, we will present the results of a local sensitivity analysis of the parameters on the threshold quantity $(R_0)$. Finally, Section 8 contains the concluding remarks.

2.   A two-stage-structure model in the human population

The proposed model is of the host–vector type, using ordinary differential equations (ODE) to describe the transmission of dengue; the population is considered to be constant.

The model has six compartments for the human population: susceptible children ($S_C$), susceptible adults ($S_A$), infected children ($I_C$), infected adults ($I_A$), recovered children ($R_C$), and recovered adults ($R_A$). The vector population has two compartments: susceptible ($S_V$) and infected ($I_V$).

The total population of humans is given by $N_h$, and $\mu_h$ is the rate of natural mortality. The parameter $\eta$ denotes the rate at which children become part of the adult population. Since dengue infection is brief duration compared to the time a child reaches adulthood, for this reason, we do not incorporate the rate of maturation $\eta$ from infected children to infected adults. We denote the contact rate that a virus-carrier mosquito has of susceptible children and adults by $\beta_{h}/N_{h}$ and the incidence of new infections transmitted by the mosquitoes are given again by $\frac{\beta_{h}S_{C}I_{V}}{N_h}$ and $\frac{\sigma\beta_{h}S_{A}I_{V}}{N_h}$, respectively.

The parameter $\sigma$ is the ratio between the probability that an adult is bitten by a mosquito and the probability that a child is bitten by a mosquito: it is the relative risk, and is assumed to be between zero and unity. Susceptible mosquitoes start carrying the virus after getting into contact with an infected human at a rate of $\beta_{v}/N_{h}$, so that the incidence of newly infected mosquitoes per infected children and adults is given by $\frac{\beta_{v}S_{V}I_{C}}{N_h}$ and $\frac{\sigma \beta_{v}S_{V}I_{A}}{N_h}$, respectively. The recovery rate for children and adults are given by $\gamma_C$ and $\gamma_A$, respectively. The total mosquito population $N_v$ is constant, and $\mu_v$ is the natural mortality rate of the mosquito. This dynamics can be represented by the compartmental diagram of Figure 1.

This model can be represented by the following system of ODEs:

where it is assumed that $S_C + S_A + I_C + I_A + R_C + R_A = N_h$ and $S_V + I_V = N_V$ are the sizes of the human and vector populations, respectively.

Since $S_V = N_V - I_V$ and substituting in the last equation of (2.1), and since the variables $R_{C}(t)$ and $R_{A}(t)$ do not appear in the equations for $S_{C}(t)$, $I_{C}(t)$, $S_{A}(t)$ and $I_{A}(t)$ of system (2.1), it is sufficient to consider the behavior of the following system:

The feasible region of system (2.1) is given by,

3.   Preliminaries

3.1.   Existence of equilibrium

Setting the equations of (2.2) equal to zero and expressing the variables $S_C, S_A, I_C$ and $I_A$ in terms of $I_V$, the equilibria of this system are obtained. Solving for $S_{C}$ and $S_{A}$ of the first and second equations, respectively, of (2.2), we obtain:

Solving for $I_C$ and $I_A$ of the third and fourth equations in (2.2), respectively, we obtain:

Assuming that $I_V = 0$, the disease-free equilibrium is

which always exists in $\Omega$.

To find the positive equilibrium in $\Omega$, it must be assumed that $I_V\neq0$. Substituting $I_C$ and $I_A$ in the fifth equation of (2.2), we perform various calculations, one obtains that $I_V$ satisfies the following quadratic equation:

where the coefficients are given by

where $R_0$ is a threshold value and is given by

Since all the parameters of the model are positive, it can be seen that the coefficient $a$ is negative. Accordingly, the quadratic equation (3.6) is concave down. In (3.11), all the parameters are positive, if $R_0 > 1$, then $p(0) > 0$, so that (3.6) has only one positive solution.

With these results, the following lemma can be asserted.

Lemma 3.1. Let $R_0$ be defined by (3.10). If $R_0 \leq 1$, then only the disease-free equilibrium $E_0$ exists in $\Omega$. If $R_0 > 1$, then there is a disease-free equilibrium in $\Omega$ and a unique endemic equilibrium $E^* = (S_{C}^*, S_{A}^*, I_{C}^*, I_{A}^*, I_{V}^*)$ in $int(\Omega)$ whose coordinates are given by equations (3.1)–(3.4) and the positive solution of (3.6).

3.2.   Basic reproductive number

We obtain the basic reproductive number $\hat{R}_{0}$ using the next generation method ^[25,26]. The Jacobian matrix $J$ of this subsystem at the disease-free equilibrium is decomposed as $J = F-V$, where $F$ is the transmission part and $V$ describe the transition terms associated with the model (2.2). These quantities are given, respectively, by

Now, the values of $F$ and $V$ are calculated because the basic reproductive number is given by $\rho(FV^{-1})$, where $\rho$ is defined as the spectral radius:

This threshold is $\hat{R}_{0}$, and is the square root of $R_0$, but both give the equivalent condition $R_0 > 1$ (or $R_0 \leq 1$) $\Leftrightarrow$ $\hat{R}_{0} > 1$ (or $\hat{R}_{0}\leq 1$). For the rest of this paper, we will use the threshold number $R_0$ as the basis for the analysis.

The expression for $\hat{R}_{0}$ represents the average of the secondary infections generated by a primary case when introduced in a whole susceptible human population. This can be seen as follows: an infected mosquito distributes bites in the human population for the rest of its life, and a proportion $\frac{\beta_{h}N_{V}}{\mu_{v}N_h}$ of these bites becomes new infections. On the other hand, the number of new infections of mosquitoes by infected children and infected adults during the infectious period is given by $\frac{\beta_{v}}{\eta + \mu_{h}}(\frac{\mu_{h}}{\gamma_{C} + \mu_{h}} +\frac{\sigma^{2}\eta}{\gamma_{A} + \mu_{h}})$, respectively. The geometric mean of these two quantities, which is equal to $\hat{R}_{0}$, gives the average number of secondary infections.

4.   Stability of disease-free equilibrium

In this section the stability properties of the disease-free equilibrium will be analyzed.

Theorem 4.1. $i)$ If $R_0 < 1$, then the disease-free equilibrium $E_0$ is local asymptotically stable. $ii)$ If $R_0 > 1$, then the disease-free equilibrium $E_0$ is unstable.

Proof. The Jacobian matrix of system $(2.2)$ is given by

The local stability of the disease-free equilibrium $E_0$ depends on the eigenvalues of the matrix:

where, clearly, the eigenvalues are $-\mu_h$, $-(\eta + \mu_h)$ and those obtained from the matrix

The characteristic polynomial of $P$ would be the following cubic equation,

where:

After some calculations, we get

Note that $A > 0$, and also that $R_0 < 1$ implies $B > 0$ and $C > 0$. Also,

which implies that $A\cdot B - C > 0$, and by the Routh–Hurwitz conditions, eigenvalues have negative real parts and consequently the disease-free equilibrium $E_0$ will be asymptotically stable. On the other hand, by Descarte's rule of signs, if $R_{0} > 1$, then $C < 0$ and the polynomial (4.4) will have a positive eigenvalue, in which case $E_{0}$ will be unstable.

Using a comparison theorem ^[27], we shall prove the global stability of the disease-free equilibrium for $R_0 < 1$. Following the idea used by Esteva et al. ^[28], we carry out the following proof.

Theorem 4.2. If $R_0 < 1$, then $E_0$ is globally asymptotically stable in $\Omega$.

Proof. The equations for $I_C$, $I_A$ and $I_V$ of model (2.2) are

Since $S_C\leq \frac{\mu_hN_h}{\eta + \mu_h}$, $S_A \leq \frac{\eta N_h}{\eta + \mu_h}$ and $I_V \geq 0$, the equations of (4.8) satisfy the following:

where $F$ and $V$ were defined above.

If $R_0 < 1$, then $\rho(FV^{-1}) < 1$, which is equivalent to $F-V$ having all its eigenvalues in the left half-plane. This means that the linear system $\bar{Z}' = (F-V)\bar{Z}(t)$, where $\bar{Z} = (I_C(t), I_A(t), I_V(t))$, is asymptotically stable for $R_0 < 1$, and consequently the solutions of this system tend to zero as $t$ approaches infinity. This means that the solutions of (2.2) corresponding to $I_C, I_A$ and $I_V$ approach zero as $t \rightarrow \infty$. Therefore, it can be concluded that the disease-free equilibrium is globally asymptotically stable for $R_0 < 1$ in the region $\Omega$.

5.   Conditions for the global stability of the endemic equilibrium

In the following, we consider the global asymptotic stability of a unique endemic equilibrium $E^*$. Motivated by the classic work of Korobeinikov ^[29,30] and Vargas-De-León ^[31,32] in host–vector models, we construct a Lyapunov function for an endemic equilibrium, using Volterra-type functions $W(t) = x(t)-1-\ln x(t)$.

Theorem 5.1. Let $R_0 > 1$. If

then the endemic equilibrium $E^*$ of (2.2) is globally asymptotically stable in $int(\Omega)$.

Proof. Define $L:int(\mathbb{R}_{+}^{5})\rightarrow \mathbb{R}_{+}$

The function $L(t)$ is defined, continuous and positive definite for all $S_{C}, S_{A}, I_{C}, I_{A}, I_{V}\geq 0$. Also, the global minimum $L(S_{C}, S_{A}, I_{C}, I_{A}, I_{V}) = 0$ occurs at $E^* = (S_{C}^*, S_{A}^*, I_{C}^*, I_{A}^*, I_{V}^*)$, and, therefore, $L$ is a Lyapunov function. Calculating its derivative along the solution of (2.2) we obtain

At an endemic equilibrium, we have

Using (5.2)–(5.6), we obtain

After several calculations, the following expression is obtained:

Rewriting the Korobeinikov-type functions as Volterra-type functions, we have

If condition (5.1) is satisfied, then $dL/dt \leq 0$. Note that $dL /dt = 0$ if and only if $S_{C} = S_{C}^*$, $S_{A} = S_{A}^*$, $I_{C} = I_{C}^*$, $I_{A} = I_{A}^*$ and $I_{V} = I_{V}^*$. Therefore the largest compact invariant set in $\{(S_{C}(t), S_{A}(t), I_{C}(t), I_{A}(t), I_{V}(t)): dL/dt = 0\}$ is the singleton $\{E^*\}$. If $R_0 > 1$ and $\frac{1}{\gamma_C+\mu_h}\geq \frac{\sigma}{\gamma_A+\mu_h}$, by the classical LaSalle invariance principle (Theorem 5.3 of ^[33]), the endemic equilibrium $E^*$ is globally asymptotically stable in $\Omega$. The theorem is proved.

Assuming that the recovery rate from dengue is the same for children and adults ^[34], then we can state the following corollary.

Corollary 1. Assume that $\gamma_A = \gamma_C$. If $R_0 > 1$, then the endemic equilibrium $E^*$ of (2.2) is globally asymptotically stable in $int(\Omega)$.

6.   Numerical simulations

In this section, we present the results of numerical simulations used to visualise the trajectories of model (2.2). The time courses of the populations of children and adults were obtained by numerical integration using the method Tsist5 ^[37] based on Runge–Kutta methods implemented by the DifferentialEquations.jl native package of the Julia ecosystem. We used a set of demographic and epidemiological parameters reported in the literature ^[21,35,36]. All the parameters involved in model (2.2) are summarized in Table 1.

In Figure 2, the numerical simulations show that the solutions converge to a disease-free equilibrium $E_0$. It can be seen that in the first weeks there is a greater number of cases of infected people, but then they decrease rapidly until they disappear. The values of the parameters are as in Table 1, and $\sigma = \frac{1}{3}$, $\beta_h = \frac{1}{2}$ and $\beta_v = \frac{1}{2}$. The initial conditions used are: $S_{C}(0) = 3\; 238\; 466$, $S_{A}(0) = 1\; 230\; 617$, $I_{C}(0) = 0$, $I_{A}(0) = 0$, $I_{V}(0) = 3000$. With these parameter values, $R_0 = 0.88$, this is the scenario where only a disease-free equilibrium exists, and it is locally and globally asymptotically stable according to Theorems 4.1 and 4.2, respectively.

In Figure 3, the numerical simulations illustrate that the solutions converge to an endemic equilibrium $E^*$. It can be seen that several outbreaks occur over a period of three thousand weeks, but each outbreak will have a lower number of cases than the previous one. The values of the parameters and the initial conditions are the same as in Figure 2, except for the parameters $\sigma = \frac{1}{3}$, $\beta_h = \frac{3}{4}$ and $\beta_v = \frac{3}{4}$. With these values, $R_0 = 1.99$. Since $\gamma_A = \gamma_C$, by Corollary 1, the endemic equilibrium is globally asymptotically stable.

7.   Local sensitivity analysis

In this section, a local sensitivity analysis is carried out in order to identify which parameters have the greatest influence on changes in the values of $R_0$. Following the idea used in ^[38], let $\lambda$ be any of the nine non-negative parameters: $\mu_h$, $N_h$, $\beta_h$, $\eta$, $\sigma$, $\gamma_C$, $\gamma_A$, $\beta_v$ and $N_v$ that determine $R_0$.

If a small perturbation $\delta \lambda$ is made to the parameter $\lambda$, a corresponding change in $R_0$ will occur as $\delta R_0$, where:

The normalized sensitivity index $\Psi_\lambda$ is the ratio of the corresponding normalized changes, and is defined as

The normalized sensitivity indices for each of the parameters are given by the following expressions:

We show that the normalized sensitivity indices of $\Psi_{N_h}, \Psi_{N_V}, \Psi_{\beta_h}, \Psi_{\beta_v}$ and $\Psi_{\mu_v}$ are constant, therefore they do not vary with respect to the value of the other parameters. On the other hand, the normalized sensitivity indices of $\Psi_{\mu_h}$, $\Psi_{\eta}$, $\Psi_\sigma$, $\Psi_{\gamma_{C}}$ and $\Psi_{\gamma_{A}}$ depend on the other parameters.

Using the parameters of the numerical simulations when $R_0 > 1$, we created Figure 4 and identified the influence of the parameters on the value of $R_0$.

The sensitivity indices shown in this figure indicate that the parameters $N_V$, $\beta_h$ and $\beta_v$ have perfect direct relationships with $R_0$. This indicates that increasing (or decreasing) by $1\%$ the value of the parameters will produce an increase (or decrease) by $1\%$ of the value of $R_0$. The index of the parameter $N_h$ and $\mu_v$ have a perfect inverse relationships with $R_0$. This indicates that if we increase the value of the parameter $N_h$ and $\mu_v$ by $1\%$, then the value of $R_0$ will decrease by $1\%$. The parameters $\mu_h$ and $\eta$ have the same magnitude of influence on $R_0$ but in the opposite direction. If we increase (or decrease) the value of the parameter $\mu_h$ by $10\%$, then the value of $R_0$ will increase (or decrease) by $4.9\%$, while if we increase the value of $\eta$ by $10\%$, then the value of $R_0$ will decrease by $4.9\%$. For the parameter $\sigma$, the normalized sensitivity index indicates that if its value is increased (or decreased) by $10\%$, then the value of $R_0$ will increase (or decrease) by $5.9\%$. For the parameters $\gamma_C$ and $\gamma_A$, an increase by $10\%$ in their values will decrease the value of $R_0$ by $7\%$ and $3\%$ respectively.

8.   Concluding remarks

In this paper, a mathematical model of the host–vector type was proposed. It considers the human population to be divided into two subpopulations: children and adults. This is in order to take into consideration that children are more likely to be bitten by mosquitoes than are adults.

A qualitative analysis of the model was carried out, obtaining two equilibria: a disease-free equilibrium $E_0$ and an endemic equilibrium $E^*$. It can be concluded that the stability of the equilibrium is completely determined by the value of $R_0$. If $R_0 < 1$, then the disease-free equilibrium is locally and globally asymptotically stable by Theorems 4.1 and 4.2, respectively. When $R_0 > 1$, then the disease-free equilibrium is unstable by Theorems 4.1 and an endemic equilibrium emerges in the interior of $\Omega$. A sufficient condition (5.1) for the global stability of the endemic equilibrium was obtained in Theorem 5.1. We can interpret this condition as follows: the average residence time $1/(\gamma_C+\mu_h)$ in the class of infected children is greater than or equal to the product of the relative risk of being bitten by a mosquito ($\sigma$) and the average residence time $1/(\gamma_A+\mu_h)$ in the class of infected adults. On the other hand, assuming that $\gamma_C = \gamma_A$, we obtained Corollary 1 without sufficient conditions on parameters for the global asymptotic stability of the endemic equilibrium.

Numerical simulations were carried out for the value of the parameters of the model (2.2) to observe the behavior of its solutions in the scenarios when $R_0 < 1$ and when $R_0 > 1$. Figure 1 shows that dengue will die out from the mosquito and human populations for any initial condition when $R_0 < 1$. Figure 2 shows that dengue remains present in the mosquito and human populations and will become endemic when $R_0 > 1$. Particularly, several outbreaks occur over a period of three thousand weeks, but each outbreak will have a lower number of cases than the previous one.

We observed that the relative risk ($\sigma$) of being bitten by a mosquito is quadratic in the basic reproductive number ($\hat{R_{0}}$), and $\sigma$ is less than unity. The above explains the lesser contribution of adults compared to children to the dynamics of dengue transmission; this can be seen in the simulations presented in Figure 2 (b).

A sensitivity analysis can provide insights into possible strategies for the control of an epidemic. The results of the local sensitivity analysis (Figure 4) should be considered together with simulated outputs of epidemiological models for the possible interventions to be carried out. In our two-stage-structure model for dengue, the transmission rates ($\beta_v$ and $\beta_h$) could be reduced by the use of anti-mosquito repellents or bed nets. Reducing the vector population ($N_V$) has been a classic control strategy and one of the measures with the greatest impact on vector-borne diseases. The vector mortality rate ($\mu_v$) is traditionally increased by chemical, physical, and biological agents. The recovery rates ($\gamma_C$ and $\gamma_A$) could be increased by developing antiviral therapies, which is still under investigation ^[39]. We note that the recovery rate of children $\gamma_C$ has a greater influence on $R_0$ than does the recovery rate of adults $\gamma_A$. The sensitivity analysis allowed us to determine that the relative risk of being bitten by a mosquito $\sigma$ has a greater influence on $R_0$ than do the recovery rates. So it is necessary to reduce the risks of children to mosquito bites when they spend time outdoors, maintaining self-protection such as the use of insect repellents. A systematic review conducted by Lima et al. ^[40] reports that the most effective control method was the integrated approach, considering the influence of eco-bio-social determinants in the virus–mosquito–human epidemiological chain, and community involvement. Several authors ^[40,41] have pointed out that is important to engage and educate the community as active agents of vector control, apart from achieving the implementation of integrated approach.

Finally, we believe that children's increased risk of being bitten by mosquitoes is a fact that should be incorporated into age-group dengue models, given that vaccination strategies are focusing on children.

Acknowledgments

Li-Martín is indebted to CONACYT for a fellowship that enabled her to pursue graduate studies for the degree of Doctorado en Ciencias Matemáticas. This work has been partially supported by the CONACYT-Ciencia de Frontera "Ecuaciones de evolución, sus estados estacionarios y su comportamiento asintótico con aplicaciones en Física y Biología" (project No. 684340).

Conflict of interest

The authors declare there is no conflict of interest.