The risk index for an SIR epidemic model and spatial spreading of the infectious disease

1.   Introduction

Chikungunya fever is a viral disease caused by an arboviral alphavirus transmitted to humans by Aedes mosquitoes, mainly Aedes aegypti and Aedes albopictus mosquitoes. These vectors are widely distributed in the Americas. In 1952, this virus was first isolated and described in humans in Tanzania; it was identified in the Americas in 2013, and in Mexico one year later ^[1].

Once an infected mosquito bites a person, the disease has an incubation period between 3 and 7 days. The symptoms are severe joint pain (arthralgia) and high fever (above 39℃); they also include occasional nausea, myalgia (muscle pain), vomiting and rash ^[2]. This joint pain is usually debilitating, affecting the quality of life of the patient, but the persistence of the symptoms has not been thoroughly studied.

In 2008, Dumont et al. proposed a host-vector model that considers the aquatic phase of the mosquito for Chikungunya disease; they also studied the local and global stability of the equilibrium points, fitted the model to epidemiological data from four cities in France and estimated the basic reproductive numbers of the outbreak ^[3]. Subsequently, Dumont et al developed a model that incorporates a combination of the early use of massive fumigation and mechanical control (such as the destruction of breeding sites) and concluded that it can be effective in stopping or containing the spread of Chikungunya infections with minimal environmental impact ^[4]. In 2012, Ruiz-Moreno et al. developed a climate-based stochastic model of mosquito population dynamics with an epidemiological model to identify temporal windows that carry epidemic risk they found that, in places with a marked seasonal variation in temperature, there was also an epidemic risk season that coincided with the period of the year in which mosquito populations survive and grow ^[5]. In 2019, González-Parra et al. analyzed an SEIRC-SEI model with intrinsic and extrinsic incubation periods for the dynamics of the transmission of Chikungunya that considers a constant human population and vectors. They studied the global stability of the disease-free and endemic equilibriums by using the second method of Lyapunov and they used bootstrapping and Markov chain Monte Carlo techniques to estimate some model parameters and the basic reproductive number for an outbreak in Colombia ^[6]. Abboubakar et al. included density-dependent rates and some control mechanisms in existing Chikungunya models. Their model presents a backward bifurcation, and the least squares method was used to estimate the basic reproductive number of an outbreak that occurred in Chad and Cameroon ^[7].

The phenomenon of relapse has been reported for Chikungunya infections ^{[8,9,10,11,12]}. Relapse is defined as the reappearance of arthralgia due to virus' persistence in the cells of the musculoskeletal tissue after a symptom-free period of at least one week ^[8], or after one month ^[12]. A cohort study based on data from a laboratory-based surveillance system in France was performed; the initial infection was confirmed via an antibody or polymerase chain reaction (PCR) testing. In this study, relapses of arthralgia were reported in 72% of patients; the mean number of relapses was four and the mean time between two relapses was 8 weeks ^[8]. On the other hand, a cross-sectional study of Acapulco in southern Mexico in December 2015 was performed; 66% of the population (3531 out of 5870 people) self-reported that they had been infected; 31.1% of those who suffered from Chikungunya (1098 out of 3531 people) reported at least one relapse at least 1 month after they recovered from the disease, 13% reported exactly one relapse, 12% reported two relapses, 4% reported three relapses and just 2% reported more than four relapses ^[12].

The relapse phenomenon has not been incorporated into the host-vector models for Chikungunya infections in the mathematical epidemiology literature. For this reason, the objectives of our study were to incorporate relapses into the standard model for Chikungunya, and to study the effect of the relapse rate on Chikungunya outbreaks. Finally, using the data collected in ^[12] that consisted of self-reported Chikungunya cases per month in Acapulco, Mexico, we have estimated several model parameters and the basic reproductive number of the infection.

In the present paper we address a Chikungunya model with relapse. We have made the following contributions:

● We extended the standard host-vector model of Chikungunya to incorporate the relapse phenomenon.

● The global stability of the disease-free equilibrium and the endemic equilibrium has been analyzed using the method of Lyapunov functions.

● We performed Bayesian estimates of the model parameters and the basic reproductive number of a Chikungunya outbreak in Acapulco, Mexico.

● Local sensitivity analysis was performed to measure the relative change of basic reproductive number for each parameter.

This paper is organized as follows: In Section 2, we will propose a Chikungunya model with relapse. In Section 3, we will compute the equilibria and basic reproductive number, and we will analyze the bifurcation that occurs at the disease-free equilibrium point as well as the global stability of the disease-free and endemic equilibrium points. In Section 4, we will fit the model to data from an outbreak in 2015 in Acapulco, Mexico to estimate model parameters and the basic reproductive number $(R_0)$ with the Bayesian approach by using the Hamiltonian Monte Carlo method. In Section 5, we will present the results of a local sensitivity analysis to describe the impact of the parameters on the $R_0$ value. Finally, Section 6 contains the concluding remarks.

2.   Mathematical model

Let $N_h$ and $N_v$ be the total number of humans (hosts) and the total number of mosquitoes (vectors), respectively. Both populations are divided into mutually exclusive compartments that are dependent on each individual's epidemiological state. Humans are classified into four compartments:

1. Susceptible humans ($S_h$): We assume that all humans are born susceptible to the virus; then, the population in this compartment increases according to a constant rate $\mu_h$ and decreases due to natural death at the same rate. Additionally, the size of the susceptible population decreases when a susceptible human becomes an infected human after an infected mosquito bite ($\frac{\beta_{h}b}{N_{h}}S_{h}I_{v}$); thus, the number of new infections depends on the number of susceptible humans who are bitten by a mosquito per time unit $b$ and the probability of an effective transmission from a vector to a human $\beta_{h}$. Consequently, the dynamics of the compartment of susceptible humans are given by the following equation:

2. Infected humans ($I_h$): The number of infected humans increases when a mosquito bite effectively transmits the virus to a susceptible human according to the mechanism described above, and also when symptoms return after a symptom-free period, that is, the relapse period $1/\delta_h$. The number of individuals in this compartment decreases when the symptomatology disappears (recovery rate $\gamma$) and due to natural death ($\mu_{h}$). The number of infected humans is determined by the following equation:

3. Asymptomatic humans ($A_h$): A fraction $p$ of infected individuals are assumed to undergo an asymptomatic period since a significant percentage of patients have relapsed after the initial infection. The number of asymptomatic humans increases at a rate $p\gamma$ and decreases due to natural death ($\mu_{h}$) and when symptoms occur again at a rate $\delta_{h}$. The reappearance of symptoms is due to the persistence of the virus in musculoskeletal tissue cells after a symptom-free period ($1/\delta_{h}$). This mechanism is suggested in ^[12]. As a result, the dynamics of the asymptomatic humans are described by

4. Recovered humans $(R_h)$: A fraction $(1-p)$ of infected individuals are assumed to undergo a full recovery with permanent immunity. The number of recovered humans increases at a rate $(1-p)\gamma$ and decreases due to natural death ($\mu_{h}$). This translates into the following equation:

Mosquitoes are classified into two compartments:

1. Susceptible vectors ($S_v$): Because vertical transmission has not been reported, we assume that all mosquitoes are born without the virus. Additionally, we assume that the birth and death rates are the same ($\mu_v$), and that consequently, the vector population is constant ($S_v+I_v = N_v$). As a susceptible mosquito can become infected only if it comes into contact with an individual who has the virus, it is supposed to become infected after biting an infected or asymptomatic human, because hosts in the second compartment still have a viral load in the cells of the musculoskeletal tissue. Therefore, analogous to transmission in humans, the rate of infection in vectors depends on the bites $b$ and the transmission probability from humans to vectors $\beta_v$; under the hypothesis that there is less transmission by the asymptomatic humans than by infected humans, we introduce the non-negative parameter $\kappa$ to moderate the fraction of transmission from asymptomatic humans to vectors. These assumptions give rise to the following equation:

2. Infected vectors $(I_v)$: The number of infected mosquitoes increases according to the transmission mechanism described above and it decreases when infected mosquitoes die, since a mosquito never recovers and its life expectancy is not modified. Therefore, the number of infected mosquitoes is subject to the following equation:

In summary, the dynamics between the hosts and vectors that determine the spread of Chikungunya (see Figure 1), including relapse, translate into the following system of ordinary differential equations:

Figure 1.  Schematic diagram of Chikungunya transmission model with relapse.

DownLoad: Full-Size Img PowerPoint

Here, both total populations are considered constant, so $N_h = S_h+I_h+A_h+R_h$ y $N_v = S_v+I_v$. Under these conditions, system (2.1) can be written as follows:

3.   Equilibria, basic reproductive number, bifurcation analysis, and global stability properties

We analyze the model given in (2.2) in the following epidemiologically feasible region:

System (2.2) contains two epidemiologically feasible equilibrium points in the non-negative orthant $\mathbb{R}_{+}^{4}$ according to direct calculation: the disease-free equilibrium $P_0 = (N_h, 0, 0, 0)$ and a unique endemic equilibrium $P^* = (S^*_{h}, I^*_{h}, A^*_{h}, I^*_{v})$, with

where

For the computation of the basic reproductive number, we use the next-generation operator method introduced in ^[13]. The notation is as follows: $F$ is a non-negative matrix of the new transmission terms, and $V$ is an $M$-matrix of the transition terms for individuals between compartments. For system (2.2), the matrices $F$ and $V$ are, respectively,

It follows that the basic reproductive number, denoted by $R_{0} = \rho(FV^{-1})$, where $\rho$ is the spectral radius, is given by

In the following subsections, without loss of generality, we present the bifurcation and global stability results in terms of the threshold quantity $\hat{R}_{0}$.

3.1.   Forward bifurcation

It is easy to prove that the disease-free equilibrium $P_{0}$ is non-hyperbolic when $\hat{R}_{0} = 1$ since one of its eigenvalues vanishes; consequently, in this case, a bifurcation could occur. This subsection shows that this is indeed the case. For this purpose, Theorem 4.1 of ^[14] will be used.

First, the abbreviated set of equations in (2.2) is rewritten in the following form:

to which the following new variables have been introduced:

If we assume that $\beta _{v} = \xi \beta _{h}$ and $\phi \equiv \beta _{h}b$ is the bifurcation parameter, system (3.4) can take the following form:

where we have made the following identifications:

and the components of the vector field of (3.5) have been denoted by $f_{i}$, with $i = 1, ..., 4$. When $\hat{R}_{0} = 1$, we write $\phi = \phi ^{\ast }$, which, according to (3.2), satisfies the following relationship:

The Jacobian matrix of system (3.5), evaluated at the equilibrium point $P_{0}$ when $\phi = \phi ^{\ast },$ is given as

It can easily be shown that this matrix has a zero eigenvalue.

According to the procedure described in ^[14], the eigenvector corresponding to this zero eigenvalue is

while the left eigenvector corresponding to the same eigenvalue of the transpose of the matrix in (3.8) is

Both vectors satisfy the condition that $w\cdot v = 1$. The nonzero second partial derivatives of the components $f_{i}$ are given by

and

Thus, the quantities $a$ and $b$ as indicated in Theorem 4.1 of ^[14] and written in terms of the previous derivatives, are, respectively,

and

where $w_{1}$, $w_{2}$, $w_{3}$ and $w_{4}$ are the first, second, third and fourth components of the eigenvector in (3.9); meanwhile, $v_{2}$ and $v_{4}$ are the first and fourth components of the left eigenvector given in (3.10).

It can be shown that, according to the Routh-Hurwitz stability criterion, to guarantee the local asymptotic stability of $P_{0}$ when $\hat{R}_{0} < 1$, the following relations must be satisfied:

and

Based on relations (3.7), (3.15) and (3.16), it can be shown that $w_{1} < 0, \; w_{2} > 0, \; w_{3} > 0$ and $w_{4} > 0$, while $v_{1} > 0, \; v_{2} > 0$ and $v_{4} > 0.$ Thus, according to (3.13) and (3.14), we have that $a < 0$ and $b > 0.$ Therefore, since $a < 0$ and $b > 0$, part iv of Theorem 4.1 in ^[14] must hold. Alternatively, this result can be formulated as follows.

Theorem 3.1. The disease-free equilibrium point $P_{0} = (N_h, 0, 0, 0)$, when $\hat{R}_{0} = 1,$ presents a forward bifurcation.

As a consequence of this result, when $\hat{R}_{0} > 1$, there is a family of asymptotically stable infected equilibrium points, which we will denote as $P^* = \left(S_{h}^{\ast }, I_{h}^{\ast }, A_{h}^{\ast }, I_{v}^{\ast }\right)$, constituting the upper branch of these types of bifurcations.

3.2.   Global stability of the disease-free equilibrium

We obtained some conditions on the global stability of the disease-free equilibrium of system (2.2) by using the method of Lyapunov functions. The usual process to construct a Lyapunov function for a disease-free equilibrium in epidemic or intra-host viral infection models is to introduce the Volterra-type function to the susceptible compartment and linear functions to the other compartments, and to determine the constants to guarantee the negativity of the derivative of the Lyapunov function along trajectories.

Theorem 3.2. If $\hat{R}_{0}\leq 1$, then the disease-free equilibrium $P_0 = (N_h, 0, 0, 0)$ of system (2.2) is globally asymptotically stable in $\Delta$.

Proof. We construct the following Lyapunov function for system (2.2):

The function $V(t)$ is defined, continuous and positive definite for all $S_{h}, I_{h}, A_{h}, I_{v}\geq 0$. Additionally, the global minimum $V(S_{h}, I_{h}, A_{h}, I_{v}) = 0$ occurs at $P_0 = (N_h, 0, 0, 0)$, and, therefore, $V$ is a Lyapunov function. By calculating its derivative along the solution of (2.2), we obtain

If $\hat{R}_0\leq1$, then $dV/dt\leq0$. Note that $dV /dt = 0$ if and only if $S_{h} = N_{h}$, $I_{h} = 0$ and $I_{v} = 0$ or if $\hat{R}_0 = 1$, $S_{h} = N_{h}$ and $I_{v} = 0$. Therefore, the largest compact invariant set in $\{(S_h, I_h, A_h, I_v): dV/dt = 0\}$ is the singleton $\{P_0\}$. By the classical LaSalle invariance principle (Theorem 5.3 of ^[15]), $P_0$ is globally asymptotically stable in $\Delta$ if $\hat{R}_0\leq1$.□

3.3.   Global stability of the endemic equilibrium point

We consider the global asymptotic stability of a unique endemic equilibrium $P^*$ by using linear combinations of Volterra-type functions. The usual process to construct a Lyapunov function for a positive equilibrium in epidemic ^{[16,17,18,19]} or intra-host viral infection ^[20] models is to propose the use of the Volterra-type function, and to determine the constants to guarantee the negativity of the derivative of the Lyapunov function along trajectories.

Theorem 3.3. If $\hat{R}_{0} > 1$, then the endemic equilibrium $P^* = \left(S_{h}^{\ast }, I_{h}^{\ast }, A_{h}^{\ast }, I_{v}^{\ast }\right)$ of system (2.2) is globally asymptotically stable in $int(\Delta)$.

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

The function $L(t)$ is defined, continuous and positive definite for all $S_h, I_h, A_h, I_v\geq 0$. Additionally, the global minimum $L(S_h, I_h, A_h, I_v) = 0$ occurs at $P^* = (S_h^*, I_h^*, A_h^*, I_v^*)$, and, therefore, $L$ is a Lyapunov function. We calculate the time derivative of $L(t)$:

The coordinates $S^*_{h}$, $I^*_{h}$, $A^*_{h}$ and $I^*_{v}$ of the endemic equilibrium (3.1) satisfy the following equations:

Using the identities given in (3.18), we have

After several calculations, we have

Therefore,

The arithmetic mean is greater than the geometric mean, i.e., the terms $\left(2-\frac{S_h^*}{S_h}-\frac{S_{h}}{S^*_h}\right)$, $\left(2-\frac{I_v^*}{I_v}-\frac{I_{v}}{I_v^*}\right)$, $\left(3-\frac{S_h^*}{S_h}-\frac{I_v^*I_{h}}{I_vI^*_{h}}-\frac{I^*_{h}S_{h}I_{v}}{I_{h}S^*_hI^*_{v}}\right)$, $\left(2-\frac{A_{h}^*I_h}{A_{h}I^*_h}-\frac{A_{h}I_h^*}{A^*_{h}I_h}\right)$ and $\left(4-\frac{S_h^*}{S_h}-\frac{I_v^*A_{h}}{I_vA^*_{h}}-\frac{A_h^*I_{h}}{A_hI_{h}^*}-\frac{I^*_{h}S_{h}I_{v}}{I_{h}S^*_hI^*_{v}}\right)$ are negative. Consequently, $dL/dt \leq 0$ for any coordinate values $(S_h, I_h, A_h, I_v)$ and $dL /dt = 0$ if and only if $S_{h} = S_{h}^*$, $I_{h} = I_{h}^*$, $A_{h} = A_{h}^*$ and $I_{v} = I_{v}^*$. Therefore, the largest compact invariant set in $\{(S_h, I_h, A_h, I_v): dL/dt = 0\}$ is the singleton $\{P^*\}$. If $\hat{R}_{0} > 1$, by the classical LaSalle invariance principle (Theorem 5.3 of ^[15]), the endemic equilibrium $P^*$ is globally asymptotically stable in $int(\Delta)$.□

4.   Bayesian approach

4.1.   Bayesian model

We consider the model

where $y_i$ is the i-th observation of the data that represents the number of infected humans, $\boldsymbol{\theta} = (\beta_h, \beta_v, b, \mu_v, \gamma, p, \delta_h, \kappa, N_v)$ is the unknown vector of parameters of the model in (2.2) to estimate, $X_{\boldsymbol{\theta}}$ is the numerical solution of the system (2.2) for the state variable $I_h$ according to the Runge–Kutta fourth-order method and $\varepsilon(t_i)$ is the random error. The Bayesian statistical model is used in situations in which the dynamics of the expected value $E(y_i) = X_{\boldsymbol{\theta}}(t_i)$ are described by the host-vector model and there is an error around it, i.e, $y_i - E(y_i) = \varepsilon(t_i)$. The errors in each measurement are assumed to be independent even if the expected values $E(y_i)$ follow the host-vector model.

We estimate $\boldsymbol{\theta}$ by using Bayesian inference, which allowed us to add previous information from the literature to the data to generate the distribution for each parameter in $\boldsymbol{\theta}$. According to Bayes' theorem,

where $P(\boldsymbol{\theta})$ is the prior distribution, $P(\boldsymbol{y}|\boldsymbol{\theta})$ is the likelihood function and $P(\boldsymbol{\theta}|\boldsymbol{y})$ is the posterior distribution. We assume that $\varepsilon(t_i)\sim N(0, \sigma^2)$ and $y_i \sim N(X_{\theta}(t_i), \sigma^2)$. The likelihood function is given by

To choose the prior distributions, we conducted a literature review of the parameters $\beta_h$, $\beta_v$, $b$, $\mu_v$ and $\gamma$ reported in host-vector models for Chikungunya ^[3,4,5,6,7], the relapse period $1/\delta_h$ from the epidemiological literature ^[8,12] and the parameters of the town (Acapulco, Mexico) where the epidemic outbreak occurred ($\mu_h$, $N_h$ and $N_v$) ^[12,21]. The parameter values are given in the third column of Table 1. We used uniform densities for most of the parameters, whose support includes the values collected from the literature. We only use beta informative prior distributions for the transmission probability parameters. The prior distributions are given in the fourth column of Table 1. Cole ^[22] showed that the lack of identifiability results in a strong dependence on the prior information, and that, if you use informative prior information, the estimates will be close to the true parameter values.

Table 1.  A review of the values of the model parameters and the selected prior distributions. For the beta distributions $Beta(\alpha, \beta)$, $\alpha$ and $\beta$ are the shape parameters. For the uniform distributions $U(a, b)$, $a$ and $b$ are the minimum and maximum values.

Parameter	Description	Reference	Mean value or range of values	Prior distribution
$\beta_h$	Transmission probability: from vector to human	[7]	0.9999 [0.6, 1]	$Beta(5, 2)$
		[4]	0.375
		[3]	[0.5, 0.8]
		[5]	0.67 [0.26, 1]
$\beta_v$	Transmission probability: from human to vector	[7]	0.6 [0.6, 1]	$Beta(5, 2)$
		[4]	0.375
		[3]	0.37
$b$	Number of bites	[7]	2.4676 [1, 3]	$U(0, 4)$
		[4]	1
		[3]	0.5 or 1
$\mu_v$	Birth and death rates of mosquitoes $[month^{-1}]$	[4]	2.72	$U(2, 4.5)$
		[3]	4.28
		[6]	2.143
$\gamma$	Recovery rate $[month^{-1}]$	[7]	[3.7, 4.5]	$U(2.5, 7.5)$
		[6]	[2, 6]
$p$	Fraction of infected becoming asymptomatic	Supposed		$U(0, 1)$
$\delta_h$	Relapse rate $[month^{-1}]$	[8]	0.5	$U(0, 1)$
		[12]	0.66
$\kappa$	Transmission fraction: from asymptomatic humans to vectors	Supposed		$Beta(2, 5)$
$N_v$	Total number of vectors	Supposed	2 to 4.5 times more than the number of humans	$U(11740, 26415)$
$\mu_h$	Birth and death rates of humans $[month^{-1}]$	[21]	$\frac{1}{75\times12}\approx0.0011$	Fixed value
$N_h$	Total number of humans	[12]	5870	Fixed value

---

 | Show Table

DownLoad: CSV

To sample from the posterior distributions of nine parameters, we use the No-U-Turn-Sampler algorithm, which is a Hamiltonian Monte Carlo method. The reader is referred to ^[23] for a description of the method. Three independent Hamiltonian Monte Carlo chains were initialized with a random initial value and run with 20,000 iterations. For each parameter, we calculated the potential scale reduction factor, commonly known as Rhat ^[24]. Rhat values less than or equal to 1.2 indicate convergence to a stationary distribution. Furthermore, we verified that the trace plots showed good mixing. The libraries DifferentialEquations ^[25] and turing ^[26] of the Julia software ^[27] were used.

The Bayesian point estimators were the means of the posterior distributions corresponding to a minimization of the expected squared error loss function. The 95% credible interval (CrI) was calculated by using the 2.5th and 97.5th percentiles.

4.2.   Bayesian estimation

In Table 2, we report the point and interval estimations of parameters and the basic reproductive number. We will discuss the convergence criteria of the Markov chain Monte Carlo method in the Appendix. The estimated fraction of infected individuals who become asymptomatic, $p$, was approximately 0.6536, and with a probability of 0.95, this value is between 0.3916 and 0.8324. The estimated fraction of transmission from asymptomatic humans to vectors was approximately 0.2558, and with a probability of 0.95, this value is between 0.0385 and 0.6030. The estimated relapse rate $\delta_h$ was 0.7506 (95% $CrI$: 0.2914–0.9906). The estimated basic reproductive number $R_0$ was 2.80 (95% $CrI$: 2.65–4.50). In Figure 2, we show the fit of the model to the number of humans infected with Chikungunya and numerical simulations of the susceptible human, infected human, asymptomatic human and infected vector compartments as performed by using the Bayesian mean estimation of the parameters.

Table 2.  Results from the Bayesian estimation of the model parameters and basic reproductive number, including their Rhat values.

Parameter	Mean	Median	95% Credible interval	Rhat
$\beta_h$	0.7379	0.7531	(0.4377, 0.9574)	1.0001
$\beta_v$	0.7570	0.7737	(0.4633, 0.9638)	1.0002
$b$	3.2158	3.2528	(2.2392, 3.9579)	1.0002
$\mu_v$	3.3510	3.3914	(2.0874, 4.4459)	1.0000
$\gamma$	4.5821	4.4888	(2.6747, 7.0150)	1.0002
$p$	0.6536	0.6710	(0.3916, 0.8324)	1.0001
$\delta_h$	0.7506	0.7910	(0.2914, 0.9906)	1.0000
$\kappa$	0.2558	0.2321	(0.0385, 0.6030)	1.0001
$N_v$	19440.1813	19580.1146	(12254.6677, 26040.5418)	1.0001
$R_0$	2.8057	1.8199	(2.6528, 4.5063)	*

---

 | Show Table

DownLoad: CSV

Figure 2.  Fitted model: (a) Fit of the model to the number of infected humans with Chikungunya by using means as point estimators; (b) Numerical simulations of the four compartments of humans based on a Bayesian mean estimation of the parameters. The solid dots represent the epidemiological data obtained from ^[12].

DownLoad: Full-Size Img PowerPoint

5.   Local sensitivity analysis

We performed a sensitivity analysis to measure the dependency of a variable on the parameters and thus identify the parameters that have the highest impact. The sensitivity index approximates the fractional change in a variable $X$ that results from a unit fractional change in a parameter $\alpha$ when the other parameters are kept constant, and it is given by $E_\alpha = \frac{\alpha}{X}\frac{\partial X}{\partial \alpha}$ ^[28]. If $E_\alpha < 0$, the relationship between the parameter $\alpha$ and $X$ is inversely proportional. Otherwise, if $E_\alpha$ is positive, the relationship is directly proportional.

Using the basic reproductive number $R_0$, as determined by using (3.3), we calculated the sensitivity index for each parameter:

We have used the fixed values of $\mu_h$ and $N_h$ from Table 1 and the Bayesian mean estimation of the model parameters from Table 2 to create Figure 3. The sensitivity indices shown in this figure indicate that the parameters $\beta_h$, $\beta_v$, $b$, $p$, $\kappa$ and $N_h$ have direct relationships with $R_0$. This implies that decreasing the values of these parameters will reduce $R_0$. On the contrary, the parameters $\mu_v$, $\gamma$, $\delta_h$, $\mu_v$ and $N_h$ have inverse relationships with $R_0$. This implies that increasing the values of these parameters will decrease $R_0$. From Figure 3, we observe that the fraction of infected individuals who become asymptomatic, $p$, has the greatest (direct or inverse) impact on $R_{0}$, followed by the number of bites $b$, the transmission probability ($\beta_h$ and $\beta_v$), the total number of vectors $N_v$ and the death rate of mosquitoes $\mu_v$.

Figure 3.  Local sensitivity analysis of the basic reproductive number $R_0$ using the Bayesian mean estimation of the parameters.

DownLoad: Full-Size Img PowerPoint

Among the new parameters incorporated into the Chikungunya model, we note that the relapse rate $\delta_h$ and the fraction of transmission from asymptomatic humans to vectors $\kappa$ have a weak impact on $R_0$. We interpret the negative relationship between $\delta_h$ and $R_0$ as increasing the relapse rate, which decreases $R_0$, that is to say, by decreasing the relapse period ($1/\delta_h$), $R_0$ is reduced. On the other hand, the fraction of infected individuals who become asymptomatic, $p$, has a strong impact on $R_0$.

6.   Concluding remarks

Mathematical models for Chikungunya have incorporated the aquatic phase of the mosquito ^[3,29], control mechanics ^[4,7] and the incorporation of intrinsic and extrinsic incubation periods ^[6,29]. The phenomenon of relapses in Chikungunya virus infections has been reported in many studies ^{[8,9,10,11,12]}, but it has not been considered in mathematical models. In this paper, we proposed the first model that incorporates the recurrence of symptoms into the dynamics of Chikungunya.

First, we carried out a qualitative analysis of the Chikungunya model with relapses. The system has the classic equilibrium points, namely, disease-free and endemic equilibrium points. We calculated the basic reproductive number ($R_{0} = \sqrt{\hat{R}_{0}}$) for the disease by using the next-generation operator method. We proved that the disease-free equilibrium point $P_{0}$, when $\hat{R}_{0} = 1,$ presents a forward bifurcation. We analyzed the global stability of the equilibrium points by using the Lyapunov function method. We proved that the disease-free equilibrium is globally asymptotically stable when $\hat{R}_{0}\leq1$, and that the endemic equilibrium is globally asymptotically stable when $\hat{R}_{0} > 1$ in the interior of the feasible region.

Second, a Chikungunya outbreak occurred in Acapulco, Mexico, and we used the epidemiological data of self-reported cases from this outbreak obtained in ^[12] to fit the state variable $I_h$ of the model. We used the Bayesian approach to estimate nine model parameters and two fixed parameter values. For this outbreak, we estimated that the fraction of infected individuals who become asymptomatic, $p$, is between 0.3916 and 0.8324. These estimates of the infected people who experience relapses with rheumatoid symptoms are consistent with those reported in the literature which ranges from approximately 30% ^[12] to 72% ^[8]. The relapse period ($1/\delta_h$) is estimated to be between 1.00 and 3.43. Infected asymptomatic humans play a relatively small role in the transmission of the infection; the fraction of transmission from asymptomatic humans to vectors ranged from 0.0385 to 0.6030. We estimated that the basic reproductive number was equal to 2.80 (95% $CrI$: 2.65–4.50); this value plays an important role as a bifurcation value and in epidemiological interpretation. Recently, Haider et al. ^[30] estimated, using a frequentist meta-analysis, the basic reproductive number for Chikungunya based on various studies that have estimated this value; they obtained a value of 3.4 (95% confidence interval: 2.4–4.2). Our Bayesian estimate of $R_0$ is very similar to that reported in ^[30].

Third, we performed a local sensitivity analysis of the basic reproductive number by using the Bayesian mean estimation of the model parameters. Among the new parameters incorporated into the Chikungunya model, we found that the fraction of infected individuals who become asymptomatic, $p$, has a strong impact on $R_0$. The fraction $p$ and the relapse period $1/\delta_h$ could be reduced with antiviral therapies, and, consequently, $R_0$ could be decreased. Conversely, the relapse rate $\delta_h$ has a weak impact on $R_0$, but it is not insignificant.

Fourth, to assess the efficiency of insecticide control in the presence and absence of relapses, we varied the mosquito death rate and examined its impact on $R_0$. The results for these two scenarios are shown in Figure 4. For the case in which the fraction $p$ is considered, we note that we need a higher mortality rate to ensure that the disease is controlled, with an $R_0$ value below unity. When relapses are not considered in the modeling, it intuitively follows that less control is needed over the mosquito population. This demonstrates the importance of incorporating relapse into Chikungunya models for the estimation of the dynamics of the disease and the study of control measures and protocols.

Figure 4.  Values of the basic reproductive number obtained by varying the death rate of mosquitoes $\mu_v$. On the left side, the fraction of infected individuals who become asymptomatic ($p$) is considered, and, on the right side, the fraction $p$ is not considered.

DownLoad: Full-Size Img PowerPoint

Finally, this work has some limitations. The occurrence of Chikungunya cases was self-reported through a designed questionnaire, and the cases were not confirmed by PCR or serological tests. Our model considers that the relapse rate is constant, but the epidemiological literature reports variability in the relapse periods that should be considered in future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was partially supported by the CONAHCYT-Ciencia de Frontera project "Ecuaciones de evolución, sus estados estacionarios y su comportamiento asintótico con aplicaciones en Física y Biología" (project no. 684340). Vázquez-Peña is indebted to CONAHCYT for the scholarship (project no. 684340) that allowed him to carry out his thesis to obtain a Bachelor's degree in Applied Mathematics at the Benemérita Universidad Autónoma de Puebla.

Conflict of interest

The authors declare that they have no conflicts of interest regarding the publication of this manuscript.

Appendix

In this Appendix, the convergence criteria of the host-vector model for Chikungunya are discussed. The chains of model parameters show good mixing, as shown in Figure 5. Table 2 shows the results of the computation of the Rhat values, thus indicating that both convergence results are consistent.

Figure 5.  On the left, the convergence chains are shown, and, on the right, the posterior distributions of the parameters of the host-vector model for Chikungunya are shown.

DownLoad: Full-Size Img PowerPoint