1. Introduction

Cholera is an acute diarrhoeal illness caused by infection of the intestine with bacterium vibrio cholerae, which lives in an aquatic organism ^[16]. The ingestion of contaminated water can cause cholera outbreaks, as John Snow proved in 1854 ^[24]. This is a possibility of transmission of the disease, but others do exist. For example, susceptible individuals can become infected if they become in contact with contaminated people. When individuals are at an increased risk of infection, then they can transmit the disease to those who live with them, by reflecting food preparation or using water storage containers ^[24]. An individual can be infected with or without symptoms and some of these can be watery diarrhea, vomiting and leg cramps. If contaminated individuals do not get treatment, then they become dehydrated, suffering acidosis and circulatory collapse. This situation can lead to death within 12 to 24 hours ^[18,24]. Some studies and experiments suggest that a recovered individual can be immune to the disease during a period of 3 to 10 years. Nevertheless, recent researches conclude that immunity can be lost after a period of weeks to months ^[19,24].

Since 1979, several mathematical models for the transmission of cholera have been proposed: see, e.g., ^{[2,3,8,9,11,14,16,17,18,19,22,24]} and references cited therein. In ^[19], the authors propose a SIR (Susceptible-Infectious-Recovered) type model and consider two classes of bacterial concentrations (hyperinfectious and less-infectious) and two classes of infectious individuals (asymptomatic and symptomatic). In ^[24], another SIR-type model is proposed that incorporates hyperinfectivity (where infectivity varies with time since the pathogen was shed) and temporary immunity, using distributed delays. The authors of ^[18] incorporate in a SIR-type model public health educational campaigns, vaccination, quarantine and treatment, as control strategies in order to curtail the disease.

Between 2007 and 2018, several cholera outbreaks occurred, namely in Angola, Haiti, Zimbabwe and Yemen ^[1,24,31]. Mathematical models have been developed and studied in order to understand the transmission dynamics of cholera epidemics, mostly focusing on the epidemic that occurred in Haiti, 2010–2011 ^[1]. In ^[16], a SIQR (Susceptible-Infectious-Quarantined-Recovered) type model, which also considers a class for bacterial concentration, is analysed, being shown that it fits well the cholera outbreak in the Department of Artibonite - Haiti, from 1$^{\text{st}}$ November 2010 until 1$^{\text{st}}$ May 2011. Furthermore, an optimal control problem is formulated and solved, where the optimal control function represents the fraction of infected individuals that are submitted to treatment through quarantine ^[16]. Many studies have been developed with the purpose to find and evaluate measures to contain the cholera spread. Nevertheless, it was not possible yet to obtain solutions in real time that can stop the cholera epidemics ^[1].

Recently, the biggest outbreak of cholera in the history of the world has occurred in Yemen ^[26]. The epidemic began in October 2016 and in February-March 2017 was in decline. However, on 27$^{\text{th}}$ April 2017 the epidemic returned. This happened ten days after Sana'a's sewer system had stopped working. Problems in infrastructures, health, water and sanitation systems in Yemen, allowed the fast spread of the disease ^[28]. Between 27$^{\text{th}}$ April 2017 and 15$^{\text{th}}$ April 2018 there were 1 090 280 suspected cases reported and 2 275 deaths due to cholera ^[32]. In ^[20], Nishiura et al. study mathematically this outbreak, trying to forecast the cholera epidemic in Yemen by explicitly addressing the reporting delay and ascertainment bias. A vaccine for cholera is currently available, but poor sanitation and the lack of access to vaccines promote the spread of the disease ^[1].

Following World Health Organization (WHO) recommendations, the first-ever oral vaccination campaign against cholera had been launched on 6$^\text{th}$ May 2018 in Yemen, but it was concluded on 15$^\text{th}$ May 2018 ^[30], due to the lack of national governmental authorization to do the vaccination ^[23]. Aid workers say that one reason of the delayed of the campaign vaccination is due to some senior Houthi officials who objected to vaccination. This campaign coincided with the rainy season and some health workers fear that this could spread the disease ^[23]. This campaign just covered four districts in Aden, which were at a high risk of fast spread of the disease, and just 350 000 individuals (including pregnant women) were vaccinated ^[23,26]. It is important to note that the vaccinated individuals represent, approximately, 1.21% of the total population, since the total population of Yemen, in 2018, is 28 915 284 ^[34]. Lorenzo Pizzoli, WHO's cholera expert, said that the campaign hoped to cover at least four million people in areas at risk (corresponding, approximately, to 14% of total population) and Michael Ryan, WHO's Assistant Director-General, revealed that they were negotiating with Yemen health authorities in order to vaccinate people from all high risks zones ^[23]. The International Coordinating Group on Vaccine Provision had planned one million cholera vaccines for Yemen in July 2017, but WHO and Yemen local authorities decided to postpone it and the doses were diverted to South Sudan. WHO and GAVI, the Vaccine Alliance, affirmed that the largest cholera vaccination is now being carried out in five countries (Kenya, Malawi, South Sudan, Uganda and Zambia). It is expected that this campaign targets more than two million people across Africa ^[23].

In this paper, we propose a SITRV (Susceptible-Infectious-Treated-Recovered-Vaccinated) type model, which includes a class of bacterial concentration. In Section 2, we formulate and explain the mathematical model. Then, in Section 3, we show that the model is mathematically well posed and it has biological meaning. We prove the existence and uniqueness of the disease-free and endemic equilibrium points and compute the basic reproduction number. The sensitivity of the basic reproduction number with respect to all parameters of the model is analysed in Subsection 3.3 and the stability analysis of equilibria is carried out in Subsection 3.4. In Section 4, we show that the model fits well the cholera outbreak in Yemen, between 27$^{\text{th}}$ April 2017 and 15$^{\text{th}}$ April 2018. Through numerical simulations, we illustrate the impact of vaccination of susceptible individuals in Yemen. We end with Section 5, by deriving some conclusions about the importance of vaccination campaigns on the control and eradication of a cholera outbreak.

We trust that our work is of great significance, because it provides a mathematical model for cholera that is deeply studied and allows to obtain important conclusions about the relevance of vaccination campaigns in cholera outbreaks. We show that, if it had existed a vaccination campaign from the beginning of the outbreak in Yemen, then the epidemic would have been extinguished. Actually, we believe that the absence of this type of prevention measures in Yemen, was one of the responsible for provoking the biggest cholera outbreak in world's history ^[26], killing 2 275 individuals until 15$^{\text{th}}$ April 2018.

---

2. Model formulation

We modify the model studied in ^[16], adding a vaccination class and considering different kinds of cholera's treatment. The model is a SITRV (Susceptible-Infectious-Treated-Recovered-Vaccinated) type model and considers a class of bacterial concentration for the dynamics of cholera. The total human population $N(t)$ is divided into five classes: susceptible $S(t)$, infectious with symptoms $I(t)$, in treatment $T(t)$, recovered $R(t)$ and vaccinated $V(t)$ at time $t$, for $t \ge 0$. Furthermore, we consider a class $B(t)$ that reflects the bacterial concentration at time $t$. We assume that there is a positive recruitment rate $\Lambda$ into the susceptible class $S(t)$ and a positive natural death rate $\mu$, for all time $t$ under study. Susceptible individuals can be vaccinated at rate $\varphi\geq 0$ and become infected with cholera at rate $\displaystyle \frac{\beta B(t)}{\kappa+B(t)}\geq 0$, which is dependent on time $t$. Note that $\beta>0$ is the ingestion rate of the bacterium through contaminated sources, $\kappa>0$ is the half saturation constant of the bacterial population, and $\displaystyle \frac{B(t)}{\kappa+B(t)}$ measures the possibility of an infected individual to have the disease with symptoms, given a contact with contaminated sources ^[18]. Any recovered and vaccinated individual can lose the immunity at rate $\omega_1\geq 0$ and $\omega_2\geq0$, respectively, becoming susceptible again. Different types of treatment for cholera infected individuals are considered based on ^[29]. The infected individuals can get a proper treatment, at rate $\delta\geq 0$, and the individuals in treatment can recover at rate $\varepsilon\geq0$. The disease-related death rates associated with the individuals that are infected and in treatment are $\alpha_1\geq0$ and $\alpha_2\geq0$, respectively. Each infected individual contributes to the increase of the bacterial concentration at rate $\eta>0$. On the other hand, the bacterial concentration can decrease at mortality rate $d>0$. These assumptions are translated into the following mathematical model:

$\begin{align} \label{ModeloColera_vaccination} &\begin{cases} S'(t) = \Lambda-\displaystyle\frac{\beta B(t)}{\kappa+B(t)}S(t) +\omega_1 R(t)+\omega_2V(t)-(\varphi+\mu)S(t), \\ I'(t) = \displaystyle\frac{\beta B(t)}{\kappa+B(t)}S(t) -\left(\delta+\alpha_1+\mu\right)I(t), \\ T'(t) = \delta I(t)-(\varepsilon+\alpha_2+\mu)T(t), \\ R'(t) = \varepsilon T(t)-(\omega_1+\mu)R(t), \\ V'(t) = \varphi S(t)-(\omega_2+\mu)V(t), \\ B'(t) = \eta I(t)-dB(t). \end{cases} \end{align}$

(2.1)

---

3. Model analysis

Throughout the paper, we assume that the initial conditions of system (2.1) are non-negative:

$\begin{equation} \label{eq:init:cond_vaccination} S(0) = S_0 \geq 0, \ I(0) = I_0 \geq 0, \ T(0) = T_0 \geq 0, \ R(0) = R_0 \geq 0, \ V(0) = V_0 \geq 0, \ B(0) = B_0 \geq 0. \end{equation}$

(3.1)

---

3.1. Positivity and boundedness of solutions

Our first lemma shows that the considered model (2.1)-(3.1) is biologically meaningful.

Lemma 1. The solutions $(S(t), I(t), T(t), R(t), V(t), B(t))$ of system (2.1) are non-negative for all $t \geq 0$ with non-negative initial conditions (3.1) in $(\mathbb{R}_0^+)^6$.

Proof. We have

$\begin{equation*} \begin{cases} \displaystyle \frac{d S(t)}{d t}\Bigg|_{\xi(S)} & = \Lambda + \omega_1 R(t) + \omega_2 V(t)> 0 \, , \\ \displaystyle \frac{d I(t)}{d t}\Bigg|_{\xi(I)} & = \displaystyle \frac{\beta B(t)}{\kappa+B(t)}S(t) \geq 0 \, , \\ \displaystyle \frac{d T(t)}{d t}\Bigg|_{\xi(T)} & = \delta I(t) \geq 0 \, , \\ \displaystyle \frac{d R(t)}{d t}\Bigg|_{\xi(R)} & = \varepsilon T(t) \geq 0 \, , \\ \displaystyle \frac{d V(t)}{d t}\Bigg|_{\xi(V)} & = \varphi S(t) \geq 0 \, , \\ \displaystyle \frac{d B(t)}{d t}\Bigg|_{\xi(B)} & = \eta I(t) \geq 0 \, , \end{cases} \end{equation*}$

where $\xi(\upsilon) = \left\{\upsilon(t) = 0 \text{ and } S, I, T, R, V, B \in C(\mathbb{R}_0^+, \mathbb{R}_0^+)\right\}$ and $\upsilon \in \{S, I, T, R, V, B\}$. Therefore, due to Lemma 2 in ^[33], any solution of system (2.1) is such that $(S(t), I(t), T(t), R(t), V(t), B(t)) \in (\mathbb{R}_0^+)^6$ for all $t \geq 0$.

Lemma 2 shows that it is enough to consider the dynamics of the flow generated by (2.1)-(3.1) in a certain region $\Omega_V$.

Lemma 2. Let

$\begin{equation} \label{eqN:1_vaccination} \Omega_H = \left\{ (S, I, T, R, V) \in \left(\mathbb{R}_0^+\right)^5 \, | \, 0 \leq S(t) + I(t) + T(t) + R(t) + V(t)\leq \frac{\Lambda}{\mu} \right\} \end{equation}$

(3.2)

and

$\begin{equation} \label{eq:maj:B} \Omega_B = \left\{ B \in \mathbb{R}_0^+ \, | \, 0 \leq B(t) \leq \frac{\Lambda\eta}{\mu d} \right\}. \end{equation}$

(3.3)

Define

$\begin{equation} \label{eq:Omega_vaccination} \Omega_V = \Omega_H \times \Omega_B. \end{equation}$

(3.4)

If $N(0) \leq \displaystyle \frac{\Lambda}{\mu}$ and $B(0) \leq \displaystyle \frac{\Lambda\eta}{\mu d}$, then the region $\Omega_V$ is positively invariant for model (2.1) with non-negative initial conditions (3.1) in $(\mathbb{R}_0^+)^6$.

Proof. Let us split system (2.1) into two parts: the human population, i.e., $S(t)$, $I(t)$, $T(t)$, $R(t)$ and $V(t)$, and the pathogen population, i.e., $B(t)$. Adding the first five equations of system (2.1) gives

$\begin{equation*} N'(t) = S'(t)+I'(t)+T'(t)+R'(t)+V'(t) = \Lambda-\mu N(t) -\alpha_1 I(t)-\alpha_2 T(t) \leq \Lambda-\mu N(t) \, . \end{equation*}$

Assuming that $N(0) \leq \displaystyle \frac{\Lambda}{\mu}$, we conclude that $N(t) \leq \displaystyle \frac{\Lambda}{\mu}$. For this reason, (3.2) defines the biologically feasible region for the human population. As it is proved in ^[16], the region (3.3) defines the biologically feasible region for the pathogen population. From (3.2) and (3.3), we know that $N(t)$ and $B(t)$ are bounded for all $t \geq 0$. Therefore, every solution of system (2.1) with initial conditions in $\Omega_V$ remains in $\Omega_V$ for all $t \geq 0$. In other words, in region $\Omega_V$ defined by (3.4), our model is epidemiologically and mathematically well posed in the sense of ^[10].

---

3.2. Equilibrium points and basic reproduction number

From now on, let us consider that $a_0 = \varphi+\mu$, $a_1 = \delta+\alpha_1+\mu$, $a_2 = \varepsilon+\alpha_2+\mu$, $a_3 = \omega_1+\mu$ and $a_4 = \omega_2+\mu$. The disease-free equilibrium (DFE) of model (2.1) is given by

$\begin{equation} \label{eq:DFE_vaccination} E^0 = \left(S^0, I^0, T^0, R^0, V^0, B^0\right) = \left(\frac{\Lambda a_4}{a_0a_4-\varphi\omega_2}, 0, 0, 0, \frac{\Lambda\varphi}{a_0a_4-\varphi\omega_2}, 0\right). \end{equation}$

(3.5)

Remark 1. Note that, because $\mu>0$, one has $a_0a_4-\varphi\omega_2 = (\varphi+\mu)(\omega_2+\mu)-\varphi\omega_2>0$.

Next, following the approach of ^[18,27], we compute the basic reproduction number $R_0$.

Proposition 3(Basic reproduction number of (2.1)). The basic reproduction number of model (2.1) is given by

$\begin{align} \label{R0_vaccination} R_0 = \frac{\beta\Lambda\eta a_4}{(a_0a_4-\varphi\omega_2)\kappa da_1}. \end{align}$

(3.6)

Proof. Consider that $\mathcal{F}_i(t)$ is the rate of appearance of new infections in the compartment associated with index $i$, $\mathcal{V}_i^+(t)$ is the rate of transfer of "individuals" into the compartment associated with index $i$ by all other means and $\mathcal{V}_i^-(t)$ is the rate of transfer of "individuals" out of compartment associated with index $i$. In this way, the matrices $\mathcal{F}(t)$, $\mathcal{V}^+(t)$ and $\mathcal{V}^-(t)$, associated with model (2.1), are given by

$\begin{align*} \mathcal{F}(t) = \left[ \begin{matrix} 0\\ \displaystyle\frac{\beta B(t)S(t)}{\kappa+B(t)}\\ 0\\ 0\\ 0\\ 0 \end{matrix} \right], \quad \mathcal{V}^+(t) = \left[ \begin{matrix} \Lambda+\omega_1 R(t)+\omega_2V(t)\\ 0\\ \delta I(t)\\ \varepsilon T(t)\\ \varphi S(t)\\ \eta I(t) \end{matrix} \right] \quad \text{ and } \quad \mathcal{V}^-(t) = \left[ \begin{matrix} \displaystyle\frac{\beta B(t)S(t)}{\kappa+B(t)}+a_0S(t)\\ a_1 I(t)\\ a_2 T(t)\\ a_3 R(t)\\ a_4 V(t)\\ d B(t) \end{matrix}\right]. \end{align*}$

Therefore, by considering $\mathcal{V}(t) = \mathcal{V}^-(t)-\mathcal{V}^+(t)$, we have that

$\begin{align*} \left[ \begin{matrix} S'(t)&I'(t)&T'(t)&R'(t)&V'(t)&B'(t) \end{matrix}\right]^T = \mathcal{F}(t)-\mathcal{V}(t). \end{align*}$

The Jacobian matrices of $\mathcal{F}(t)$ and of $\mathcal{V}(t)$ are, respectively, given by

$\begin{align*} F = \left[\begin{matrix} 0&0&0&0&0&0 \\ \displaystyle\frac{\beta B(t)}{\kappa+B(t)} & 0&0&0&0 & \displaystyle\frac{\beta\kappa S(t)}{(\kappa+B(t))^2} \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \end{matrix}\right] \text{ and } V = \left[\begin{matrix} \displaystyle\frac{\beta B(t)}{\kappa+B(t)}+a_0 & 0&0&-\omega_1&-\omega_2& \displaystyle\frac{\beta\kappa S(t)}{(\kappa+B(t))^2} \\ 0&a_1&0&0&0&0 \\ 0&-\delta&a_2&0&0&0 \\ 0&0&-\varepsilon&a_3&0&0 \\ -\varphi&0&0&0&a_4&0 \\ 0&-\eta&0&0&0&d \\ \end{matrix}\right]. \end{align*}$

In the disease-free equilibrium $E^0$ defined by (3.5), we obtain the matrices $F_0$ and $V_0$ given by

$\begin{align*} F_0 = \left[\begin{matrix} 0&0&0&0&0&0 \\ 0&0&0&0&0 & \displaystyle\frac{\beta\Lambda a_4}{(a_0a_4-\varphi\omega_2)\kappa} \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \end{matrix}\right] \quad \text{ and } \quad V_0 = \left[\begin{matrix} a_0&0&0&-\omega_1&-\omega_2& \displaystyle\frac{\beta\Lambda a_4}{(a_0a_4-\varphi\omega_2)\kappa} \\ 0&a_1&0&0&0&0 \\ 0&-\delta&a_2&0&0&0 \\ 0&0&-\varepsilon&a_3&0&0 \\ -\varphi&0&0&0&a_4&0 \\ 0&-\eta&0&0&0&d \\ \end{matrix}\right]. \end{align*}$

The basic reproduction number of model (2.1) is then given by

$\begin{align*} R_0 = \rho(F_0V_0^{-1}) = \frac{\beta\Lambda\eta a_4}{(a_0a_4-\varphi\omega_2)\kappa da_1}, \end{align*}$

found with the help of the computer algebra system $\textsf{Maple}$. This concludes the proof.

Now we prove the existence of an endemic equilibrium when $R_0$ given by (3.6) is greater than one.

Proposition 4 (Endemic equilibrium). If the basic reproduction number (3.6) is such that $R_0>1$, then the model (2.1) has an endemic equilibrium given by

$\begin{equation} \label{EndemicEquilibrium_vaccination} E^* = (S^*, I^*, T^*, R^*, V^*, B^*), \end{equation}$

(3.7)

where

$\begin{equation*} \begin{cases} S^* = \displaystyle\frac{a_1a_4\{\kappa d(a_1a_2a_3-\delta\varepsilon\omega_1) +\Lambda\eta a_2a_3\}}{\eta\tilde{D}}, \\[0.3cm] I^* = \displaystyle\frac{a_2a_3\{\beta\Lambda\eta a_4 -(a_0a_4-\varphi\omega_2)\kappa da_1\}}{\eta\tilde{D}}, \\[0.3cm] T^* = \displaystyle\frac{a_3\delta\{\beta\Lambda\eta a_4 -(a_0a_4-\varphi\omega_2)\kappa da_1\}}{\eta\tilde{D}}, \\[0.3cm] R^* = \displaystyle\frac{\delta\varepsilon\{\beta\Lambda\eta a_4 -(a_0a_4-\varphi\omega_2)\kappa da_1\}}{\eta\tilde{D}}, \\[0.3cm] V^* = \displaystyle\frac{a_1\varphi\{\kappa d(a_1a_2a_3 -\delta\varepsilon\omega_1)+\Lambda\eta a_2a_3\}}{\eta\tilde{D}}, \\[0.3cm] B^* = \displaystyle\frac{a_2a_3\{\beta\Lambda\eta a_4 -(a_0a_4-\varphi\omega_2)\kappa da_1\}}{d\tilde{D}} \end{cases} \end{equation*}$

and $\tilde{D} = a_1a_2a_3(a_0a_4-\varphi\omega_2) +\beta a_4(a_1a_2a_3-\delta\varepsilon\omega_1)$.

Proof. We note that

1. $a_1 = \delta+\alpha_1+\mu>0$, because $\delta$, $\alpha_1\geq0$ and $\mu>0$;

2. $a_2 = \varepsilon+\alpha_2+\mu>0$, because $\varepsilon$, $\alpha_2\geq0$ and $\mu>0$;

3. $a_3 = \omega_1+\mu>0$, because $\omega_1\geq0$ and $\mu>0$;

4. $\beta$, $\kappa$, $d>0$ and $\delta, \varepsilon, \varphi\geq0$;

5. $a_0a_4-\varphi\omega_2>0$ (see Remark 1);

6. $a_1a_2a_3-\delta\varepsilon\omega_1 = (\delta+\alpha_1 +\mu)(\varepsilon+\alpha_2+\mu)(\omega_1+\mu)-\delta\varepsilon\omega_1>0$, because $\alpha_1$, $\alpha_2\geq0$ and $\mu>0$;

7. $\Lambda\eta a_2a_3>0$, because $\Lambda$, $\eta$, $a_2$, $a_3>0$.

With the above inequalities, we conclude that $\tilde{D}>0$ and, consequently, that $S^*>0$ and $V^*\geq0$. The basic reproduction number is given by $\displaystyle\frac{\beta\Lambda\eta a_4}{(a_0 a_4-\varphi\omega_2)\kappa da_1}$. Thus, it follows that

$\begin{align*} &\ \beta\Lambda\eta a_4 = R_0(a_0a_4-\varphi\omega_2)\kappa da_1\\ \Leftrightarrow &\ \beta\Lambda\eta a_4-(a_0a_4-\varphi\omega_2)\kappa da_1 = R_0(a_0a_4-\varphi\omega_2)\kappa da_1-(a_0a_4-\varphi\omega_2)\kappa da_1\\ \Leftrightarrow &\ \beta\Lambda\eta a_4-(a_0a_4-\varphi\omega_2)\kappa da_1 = (a_0a_4-\varphi\omega_2)\kappa da_1(R_0-1). \end{align*}$

Therefore, we have that

$\begin{equation*} \begin{cases} I^* = \displaystyle\frac{a_1a_2a_3\kappa d(a_0a_4-\varphi\omega_2)(R_0-1)}{\eta\tilde{D}}, \\[0.3cm] T^* = \displaystyle\frac{a_1a_3\kappa d\delta(a_0 a_4 -\varphi\omega_2)(R_0-1)}{\eta\tilde{D}}, \\[0.3cm] R^* = \displaystyle\frac{a_1\kappa d\delta\varepsilon(a_0 a_4-\varphi\omega_2)(R_0-1)}{\eta\tilde{D}}, \\[0.3cm] B^* = \displaystyle\frac{a_1a_2a_3\kappa(a_0 a_4-\varphi\omega_2)(R_0-1)}{\tilde{D}}. \end{cases} \end{equation*}$

In order to obtain an endemic equilibrium, we have to ensure that $I^*$, $B^*>0$. Thus, we obtain $I^*$, $B^*>0$ if and only if $R_0-1>0\Leftrightarrow R_0>1$. In this case ($R_0>1$), we also have that $T^*$, $R^*\geq0$.

---

3.3. Sensitivity of the basic reproduction number

In this section, we are going to study the sensitivity of $R_0$ with respect to all parameters $p$ of model (2.1), computing the respective normalized forward sensitive indexes $\Upsilon^{R_0}_p$ given in Definition 1. They are presented in Table 1.

Definition 1 (See ^[6,15,25]). The normalized forward sensitivity index of a variable $z$ that depends differentiably on a parameter $p$ is defined by

$\begin{equation*} \Upsilon^z_p = \frac{\partial z}{\partial p}\times\frac{p}{|z|}. \end{equation*}$

Remark 2. When a parameter $p$ is one of the most sensitive parameters with respect to a variable $z$, then we have $\Upsilon^z_p = \pm1$. If $\Upsilon^z_p = 1$, then an increase (decrease) of $p$ by $\gamma\%$ provokes an increase (decrease) of $z$ by $\gamma\%$. On the other hand, if $\Upsilon^z_p = -1$, then an increase (decrease) of $p$ by $\gamma\%$ provokes a decrease (increase) of $z$ by $\gamma\%$ (see ^[25]).

---

3.4. Stability analysis

Now we prove the local stability of the disease-free equilibrium $E^0$.

Theorem 5 (Stability of the DFE (3.5)). The disease-free equilibrium $E^0$ of model (2.1) is

1. locally asymptotic stable, if $R_0<1$;

2. unstable, if $R_0>1$.

Moreover, if $R_0 = 1$, then a critical case occurs.

Proof. The characteristic polynomial associated with the linearised system of model (2.1) is given by

$\begin{align*} p_V(\chi) = \det(F_0-V_0-\chi I_6). \end{align*}$

In order to compute the roots of the polynomial $p_V$, we have that

$\begin{equation*} \left|\begin{matrix} -a_0-\chi&0&0&\omega_1&\omega_2& -\displaystyle\frac{\beta \Lambda a_4}{(a_0a_4-\varphi\omega_2)\kappa}\\[0.25cm] 0&-a_1-\chi&0&0&0&\displaystyle \frac{\beta \Lambda a_4}{(a_0a_4-\varphi\omega_2)\kappa} \\[0.25cm] 0&\delta&-a_2-\chi&0&0&0\\[0.25cm] 0&0&\varepsilon&-a_3-\chi&0&0\\[0.25cm] \varphi&0&0&0&-a_4-\chi&0\\[0.25cm] 0&\eta&0&0&0&-d-\chi \end{matrix}\right| = 0, \end{equation*}$

that is,

$\chi^2+(a_0+a_4)\chi+(a_0a_4-\varphi\omega_2) = 0 \vee \chi^2 +(a_1+d)\chi+a_1d-\frac{\beta\Lambda\eta a_4}{(a_0a_4-\varphi\omega_2)\kappa} = 0 \vee \chi = -a_2 \vee \chi = -a_3.$

As the coefficients of polynomial $\chi^2+(a_0+a_4)\chi+(a_0a_4-\varphi\omega_2)$ have the same sign (see Remark 1), then it follows from Routh's criterion that their roots have negative real part (see, e.g., pp.~55-56 of ^[21]). Furthermore, using similar arguments, the roots of the polynomial

$\chi^2+(a_1+d)\chi+a_1d -\frac{\beta\Lambda\eta a_4}{(a_0a_4-\varphi\omega_2)\kappa}$

have negative real part if and only if

$a_1d-\frac{\beta\Lambda\eta a_4}{(a_0a_4-\varphi\omega_2)\kappa} >0 \Leftrightarrow R_0<1.$

Therefore, the DFE $E^0$ is locally asymptotic stable if $R_0<1$; unstable if $R_0>1$; critical if $R_0 = 1$.

We end this section by proving the local stability of the endemic equilibrium $E^*$. Our proof is based on the Center Manifold Theory ^[4], as described in Theorem 4.1 of ^[5].

Theorem 6 (Local asymptotic stability of the endemic equilibrium (3.7)). The endemic equilibrium $E^*$ of model (2.1) (see Proposition 4) is locally asymptotic stable for $R_0>1$.

Proof. In order to apply the method described in Theorem 4.1 of ^[5], we are going to do the following change of variables. Let us consider that

$\begin{equation*} X = (x_1, x_2, x_3, x_4, x_5, x_6) = (S, I, T, R, V, B). \end{equation*}$

Consequently, we have that the total number of individuals is given by $N = \sum_{i = 1}^{5}x_i$. Thus, the model (2.1) can be written as follows:

$\begin{align} \label{ModeloColeraX_vaccination} \begin{cases} x_1'(t) = f_1 = \Lambda-\displaystyle\frac{\beta x_6(t)}{\kappa+x_6(t)}x_1(t) +\omega_1 x_4(t)+\omega_2x_5(t)-a_0x_1(t), \\ x_2'(t) = f_2 = \displaystyle\frac{\beta x_6(t)}{\kappa+x_6(t)}x_1(t)-a_1x_2(t), \\ x_3'(t) = f_3 = \delta x_2(t)-a_2x_3(t), \\ x_4'(t) = f_4 = \varepsilon x_3(t)-a_3x_4(t), \\ x_5'(t) = f_5 = \varphi x_1(t)-a_4x_5(t), \\ x_6'(t) = f_6 = \eta x_2(t)-dx_6(t). \end{cases} \end{align}$

(3.8)

Choosing $\beta^*$ as bifurcation parameter and solving for $\beta$, from $R_0 = 1$ we have that

$\begin{align*} \beta^* = \frac{(a_0a_4-\varphi\omega_2)\kappa da_1}{\Lambda\eta a_4}. \end{align*}$

Considering $\beta = \beta^*$, the Jacobian of the system (3.8) evaluated at $E^0$ is given by

$\begin{align*} J_0^* = \left[ \begin{matrix} -a_0&0&0&\omega_1&\omega_2 & -\displaystyle\frac{a_1d}{\eta} \\ 0&-a_1&0&0&0&\displaystyle\frac{a_1d}{\eta} \\ 0&\delta&-a_2&0&0&0\\ 0&0&\varepsilon&-a_3&0&0\\ \varphi&0&0&0&-a_4&0\\ 0&\eta&0&0&0&-d \\ \end{matrix} \right]. \end{align*}$

The eigenvalues of $J_0^*$ are obtained solving the equation $\det(J_0^*-\chi I_6) = 0$. Thus, we have that

$\begin{align*} &\ \det(J_0^*-\chi I_6) = 0\\ \Leftrightarrow&\ \chi = 0\ \vee\ \chi = -a_1-d\ \vee\ \chi = -a_2\ \vee\ \chi = -a_3\ \vee\ \chi = -\frac{1}{2}\left(a_0+a_4 \pm\sqrt{(a_0-a_4)^2+4\varphi\omega_2}\right). \end{align*}$

Note that the eigenvalue $\chi = -\frac{1}{2}\left(a_0+a_4-\sqrt{(a_0-a_4)^2 +4\varphi\omega_2}\right)$ is a negative real number, because

$\begin{eqnarray*} -\frac{1}{2}\left(a_0+a_4-\sqrt{(a_0-a_4)^2+4\varphi\omega_2}\right) & = &-\frac{1}{2}\left(\varphi+\mu+\omega_2+\mu-\sqrt{(\varphi+\mu-\omega_2-\mu)^2 +4\varphi\omega_2}\right)\\ & = &-\frac{1}{2}\left(\varphi+\omega_2+2\mu-\sqrt{\varphi^2-2\varphi\omega_2 +\omega_2^2+4\varphi\omega_2}\right)\\ & = &-\frac{1}{2}\left(\varphi+\omega_2+2\mu-\sqrt{(\varphi+\omega_2)^2}\right)\\ &\underset{\varphi+\omega_2\geq0}{ = }& -\frac{1}{2}\left(\varphi+\omega_2+2\mu-(\varphi+\omega_2)\right)\\ & = &-\mu<0. \end{eqnarray*}$

Therefore, we can conclude that a simple eigenvalue of $J_0^*$ is zero, while all other eigenvalues of $J_0^*$ have negative real part. So, the Center Manifold Theory ^[4] can be applied to study the dynamics of (3.8) near $\beta = \beta^*$. Theorem 4.1 in ^[5] is used to show the local asymptotic stability of the endemic equilibrium point of (3.8), for $\beta$ near $\beta^*$. The Jacobian $J_0^*$ has, respectively, a right eigenvector and a left eigenvector (associated with the zero eigenvalue),

$w = \left[ \begin{matrix} w_1\ \ &w_2\ \ &w_3\ \ &w_4\ \ &w_5\ \ &w_6 \end{matrix} \right]^T \quad \text{ and } \quad v = \left[ \begin{matrix} v_1\ \ &v_2\ \ &v_3\ \ &v_4\ \ &v_5\ \ &v_6 \end{matrix} \right],$

given by

$\begin{align*} w& = \left[ \begin{matrix} \displaystyle-\frac{a_4(a_1a_2a_3-\delta\varepsilon\omega_1)}{a_2 a_3(a_0a_4-\varphi\omega_2)}\ \ &1\ \ &\displaystyle\frac{\delta}{a_2} \ \ &\displaystyle\frac{\delta\varepsilon}{a_2a_3}\ \ \ & \displaystyle\frac{\varphi}{a_4}\ \ &\displaystyle\frac{\eta}{d} \end{matrix} \right]^T w_2 \end{align*}$

and

$\begin{align*} v = \left[\begin{matrix}0\ \ &1\ \ &0\ \ &0\ \ &0\ \ \ & \displaystyle\frac{a_1}{\eta}\end{matrix}\right] v_2. \end{align*}$

Remember that $f_l$ represents the right-hand side of the $l$th equation of the system (3.8) and $x_l$ is the state variable whose derivative is given by the $l$th equation, $l = 1, \ldots, 6$. The local stability near the bifurcation point $\beta = \beta^*$ is determined by the signs of two associated constants $a$ and $b$ defined by

$\begin{align*} a = \sum\limits_{i, j, k = 1}^{6}w_i w_j v_k\left[\frac{\partial^2f_k}{\partial x_i \partial x_j}(E^0)\right]_{\beta = \beta^*} \text{ and } b = \sum\limits_{i, k = 1}^{6}w_i v_k\left[\frac{\partial^2f_k}{\partial x_i \partial \phi}(E^0)\right]_{\beta = \beta^*} \end{align*}$

with $\phi = \beta-\beta^*$. As $v_1 = v_3 = v_4 = v_5 = 0$, we only have to consider the following non-zero partial derivatives at the disease free equilibrium $E^0$:

$\begin{align*} \left[\frac{\partial^2f_2}{\partial x_1\partial x_6}(E^0)\right]_{\beta = \beta^*} = \left[\frac{\partial^2f_2}{\partial x_6\partial x_1}(E^0)\right]_{\beta = \beta^*} = \frac{\beta^*}{\kappa} \ \text{ and }\ \left[\frac{\partial^2f_2}{\partial x_6^2}(E^0)\right]_{\beta = \beta^*} = -\frac{2\beta^*\Lambda a_4}{a_0a_4-\varphi\omega_2}. \end{align*}$

Therefore, the constant $a$ is

$\begin{align*} a = -\frac{2\beta^*\eta a_4}{d(a_0a_4-\varphi\omega_2)}\left( \frac{a_1a_2a_3-\delta\varepsilon\omega_1}{a_2a_3\kappa} +\frac{\Lambda\eta}{d}\right)v_2w_2^2<0. \end{align*}$

Furthermore, we have that

$\begin{equation*} \begin{split} b = v_2w_6\left[\frac{\partial^2f}{\partial x_6\partial\phi}(E^0)\right]_{\beta = \beta^*} = \frac{\Lambda\eta a_4}{\kappa d (a_0a_4-\varphi\omega_2)}v_2w_2>0. \end{split} \end{equation*}$

Thus, as

$\begin{equation*} \begin{cases} a<0\\ b>0\\ \phi = \beta-\beta^* = \displaystyle\frac{a_1\kappa d(a_0a_4 -\varphi\omega_2)}{\Lambda\eta a_4}(R_0-1)>0 \end{cases} \Leftrightarrow \begin{cases} a<0\\ b>0\\ R_0>1, \end{cases} \end{equation*}$

we conclude from Theorem 4.1 in ^[5] that the endemic equilibrium $E^*$ of (2.1) is locally asymptotic stable for a value of the basic reproduction number such that $R_0>1$.

---

4. Numerical simulations

In this section, we simulate the worst cholera outbreak that ever occurred in human history. It occurred in Yemen, from 27$^{\text{th}}$ April 2017 to 15$^{\text{th}}$ April 2018 ^[31]. As the first-ever oral cholera vaccination campaign had been launched only on 6$^{\text{th}}$ May 2018 and was concluded on 15$^{\text{th}}$ May 2018 ^[30], to describe such reality of Yemen, a numerical simulation of our model is carried out with $\varphi = \omega_2 = V(0) = 0$, that is, in absence of vaccination, and with all the other values as in Table 2. We also simulate an hypothetical situation that includes vaccination from the beginning of the outbreak, considering in that case all values of Table 2. Let us denote the numerical simulation without and with vaccination by $(NS)$ and $(NS_V)$, respectively. The curves of infected individuals for $(NS)$ and $(NS_V)$ can be observed in Figure 1, respectively in blue solid line and in blue dashed line. Our results allow us to state that if a vaccination campaign had been considered earlier in time, the number of infected individuals would have been significantly lower. Furthermore, the basic reproduction number of the simulation without vaccination is $R_0\simeq6.132305>1$ and the one with vaccination is $R_0\simeq0.753969<1$. This means that if vaccination had been considered from the beginning of the outbreak, then the spread of cholera would have been extinguished. Consequently, there would not have been so many deaths. Note that the decrease of $R_0$ with the introduction of a vaccination campaign is expected, because

$\Upsilon^{R_0}_\varphi = -\displaystyle \frac{\varphi}{\varphi+\omega_2+\mu}\simeq-0.877050<0.$

Furthermore, for $(NS)$, we obtain an endemic equilibrium point given by

$\begin{multline*} (S^*, I^*, T^*, R^*, V^*, B^*)\\ = (2.943350\times10^7, \ 1.035599 \times10^5, \ 5.954131\times10^5, \ 1.070992\times10^8, \ 0, \ 3.138180\times10^6) \end{multline*}$

and for $(NS_V)$ we have a disease-free equilibrium point given by

$(S^0, I^0, T^0, R^0, V^0, B^0) = (1.689119\times10^7, \ 0, \ 0, \ 0, \ 1.204910\times10^8, \ 0).$

Note that the previous figures correspond to the equilibrium points for the parameter values of Table 2, which can be obtained numerically for a final time of approximately 1 370 years. We also call attention to the fact that the recruitment rate $\Lambda$ of Yemen is big and this leads to a huge growth of the population.

---

5. Conclusion

In this paper, we proposed and analysed, analytically and numerically, a SITRVB model for cholera transmission dynamics. In order to fit the biggest cholera outbreak worldwide, which has occurred very recently in Yemen, we simulated the outbreak of Yemen without vaccination. Indeed, vaccination did not exist in Yemen from 27$^{\text{th}}$ April 2017 to 15$^{\text{th}}$ April 2018. Simulations of our mathematical model, with and without vaccination, show that the introduction of vaccination from the beginning of the epidemic could have changed the situation in Yemen substantially, to the case $R_0<1$, where the disease extinguishes naturally. Therefore, our research motivates and fortify the importance of vaccination in cholera epidemics.

---

Acknowledgments

This research was supported by the Portuguese Foundation for Science and Technology (FCT) within projects UID/MAT/04106/2013 (CIDMA), and PTDC/EEI-AUT/2933/2014 (TOCCATA), funded by Project 3599 - Promover a Produção Científica e Desenvolvimento Tecnológico e a Constituição de Redes Temáticas, and FEDER funds through COMPETE 2020, Programa Operacional Competitividade e Internacionalização (POCI). Lemos-Paião is also supported by the FCT Ph.D. fellowship PD/BD/114184/2016, Silva by the postdoctoral grant SFRH/BPD/72061/2010.

The authors are very grateful to an anonymous referee for reading their paper carefully and for several constructive remarks, questions and suggestions.

---

Conflict of interest

The authors declare that there is no conflicts of interest in this paper.

---