Mathematical Study for Chikungunya Virus with Nonlinear General Incidence Rate

Abstract

In this article, we examine the dynamics of a Chikungunya virus (CHIKV) infection model with two routes of infection. The model uses four categories, namely, uninfected cells, infected cells, the CHIKV virus, and antibodies. The equilibrium points of the model, which consist of the free point for the CHIKV and CHIKV endemic point, are first analytically determined. Next, the local stability of the equilibrium points is studied, based on the basic reproduction number (${\mathcal{R}}_0$) obtained by the next-generation matrix. From the analysis, it is found that the disease-free point is locally asymptotically stable if ${\mathcal{R}}_0\le 1$, and the CHIKV endemic point is locally asymptotically stable if ${\mathcal{R}}_0>1$. Using the Lyapunov method, the global stability analysis of the steady-states confirms the local stability results. We then describe our design of an optimal recruitment strategy to minimize the number of infected cells, as well as a nonlinear optimal control problem. Some numerical simulations are provided to visualize the analytical results obtained.

1. Introduction

The Chikungunya virus (CHIKV) belongs to the Togaviridae family of the genus Alphavirus [1,2]. It takes the form of an enveloped spherical viral particle with a diameter of 65 nm, containing a single RNA of positive polarity, encoding two polyproteins. The RNA is directly infectious and, therefore, serves as both a genome and messenger RNA. The 11,000 to 12,000 bases of its genome ultimately allow the synthesis of nine proteins, obtained after cleavage of the polyproteins by viral and cellular proteases. At the end of the viral cycle, which includes protein synthesis, replication, and assembly of viral particles, the viral proteins bud on the plasma membrane of the infected cell, and then penetrate the membrane. The virus is transmitted to humans when anticoagulant saliva is injected into the blood of the person bitten by an infected mosquito [1,2]. We now know that the virus does not infect circulating blood cells, but rather, macrophages and adherent cells (endothelial, epithelial, fibroblast).

In [3], Sourisseau et al. first adapted tools (flow cytometry, immunofluorescence, electron microscopy, etc.) to specify and quantify the CHIKV. Therefore, they proved in vitro that CHIKV does not replicate in circulating blood cells (lymphocytes, monocytes), but that it replicates in macrophages (phagocytic cells which are of blood origin but localized in tissues). These cells itself infect tissues, such as muscles and joints. The CHIKV also infects the so-called “adherent” cells, such as endothelial cells, epithelial cells, and fibroblasts. A study conducted by Ozden at al. [4] has shown that, in infected people, certain cells present in muscle tissue are targets of the CHIKV. Their work was based on the study of patient biopsies. They found that, in a biopsy taken at an acute stage of the disease from one patient, and in another taken at a later stage from another patient, the precursor cells of muscle cells—the satellite cells—were infected with the virus. In addition, these cells have been found, in cell culture, to be very permissive to the virus. The authors are now trying to discover if these cells play the role of a “reservoir” of the virus, which would explain the recurrence of muscle pain observed in some patients. [5] stated that like any virus, the CHIKV enters its host’s cells and uses the cellular machinery there to replicate. In particular, it is shown that this arbovirus replicates in macrophages, cells of the immune system specializing in phagocytosis and present in several tissues and organs. On the other hand, it does not reproduce in the T and B lymphocytes which are circulating in the blood. CHIKV also infects epithelial cells (organ wall), endothelial cells (inner wall of blood vessels), as well as cells of fibrous tissues (fibroblasts).

The contribution of mathematics is made initially through modelling [6,7,8,9,10,11,12,13,14]. This stage of mathematical modelling makes it possible to test, without wasting time or expense, the control measures that are envisaged: preventive measures, isolation of patients, treatments, vaccinations, and so forth. The model, nevertheless, does not completely reflect reality and is not intended to reproduce it in full. It must reproduce the characteristics of the phenomenon studied according to the objectives set for the framework of the study as well as possible. Modelling such a phenomenon consists of applying mathematical tools to a fragment of reality. It is transforming a need into equations, trying as much as possible to account for the constraints identified. Modelling is the most delicate, the longest, and often the most perilous step. Indeed, it is necessary to successfully understand the real problem to try to propose a suitable model. The first proposed attempt only very rarely meets expectations, and several modifications then follow, until a model is reached that groups together and reflects the maximum number of constraints that the real phenomenon must observe. If this step is neglected or omitted, if the constraints are not well-posed, then one ends up with a mathematical formulation which does not correspond to the problem. The resolution of the mathematical problem then provides a solution not suited to the concrete problem. However, if the problem is well-posed, then the next step is to solve the problem, that is, to analyse the model to understand, predict, and act.

Researchers have proposed several mathematical systems describing the CHIKV dynamics (see, e.g., [15,16,17,18,19,20,21]). Most of the proposed models were designed to describe the disease transmission from mosquitoes to human populations. In [22], Wang and Liu proposed a mathematical model which assumed that the contamination of a monocyte (a type of leukocyte, or white blood cell) occurs when a monocyte comes into contact with the virus. Several mathematical models have been proposed using both cellular and viral infections [13,20,23,24,25]. Elaiw et al. in [17] proposed and studied a mathematical model with two routes of infection. Later, Elaiw et al. [18,19] studied both routes of infection in a CHIKV dynamics model with Holling type-II incidence rates. Several forms of incidence rates have been applied in epidemic models [26].

In the present paper, we reconsider the mathematical system given in [18] by Elaiw et al. by considering two main changes relevant to an applied perspective:

1.

Considering the adherent cells as the main target for CHIKV and not monocytes, as mentioned in a series of previous papers [15,16,17,18,19,20].

2.

Considering general nonlinear increasing incidence rates with respect to the uninfected cells. This choice is motivated by the fact that the number of effective contacts between infective cells and susceptible cells or between the virus and susceptible cells may increase at high infective levels due to crowding of susceptible cells.

The basic reproduction number, ${\mathcal{R}}_0$, was determined based on the next-generation matrix method [27,28,29]. Local stability analysis was carried out using local linearisation. Global stability analysis of the steady-states was studied using Lyapunov stability. We proved that the CHIKV-free equilibrium point, $E_0$, is globally asymptotically stable if ${\mathcal{R}}_0\le 1$ and the infected equilibrium point, $E_1$, is globally asymptotically stable if ${\mathcal{R}}_0>1$. Furthermore, we proposed an optimal recruitment strategy to optimize the number of infected cells. We designed an optimal control problem. The resolution of the state system was reached by applying an improvement to the Gauss-Seidel-like implicit finite-difference method. The adjoint system was solved using a first-order backward-difference. Some numerical simulations confirming the obtained results are given.

2. Mathematical Model and Results

Consider the following mathematical model describing the CHIKV dynamics with two routes of infection which is a generalization of the mathematical model given in [18].

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\dot{s}\hfill & =\mathrm{\Gamma }-\delta s-{\eta }_{1}pf\left(s\right)-{\eta }_{2}yf\left(s\right),\hfill \\ \dot{y}\hfill & ={\eta }_{1}pf\left(s\right)+{\eta }_{2}yf\left(s\right)-ϵy,\hfill \\ \dot{p}\hfill & =\pi y-cp-rxp,\hfill \\ \dot{x}\hfill & =\lambda +\rho xp-mx.\hfill \end{array}\right.\end{array}$$

(1)

The variables $s,y,p$ and x describe the number of uninfected and infected cells, the CHIKV virus, and antibodies, respectively. The model parameters are nonnegative and are described as follows.

Parameter	Description
$\mathrm{\Gamma }$	uninfected cells recruitment rate
$\delta$	Uninfected cells mortality rate
$ϵ$	Infected cells mortality rate
c	CHIKV mortality rate
m	Antibodies loss rate
$\pi$	Generation rate of the virus by infected cells
r	Rate at which Antibodies attack the virus
$\lambda$	Antibodies expansion rate
$\rho$	Proliferate rate of antibodies, once antigen is encountered
${\eta }_1,{\eta }_2$	Incidence rate constants

Note that $({\eta }_{1}pf\left(s\right)+{\eta }_{2}yf\left(s\right))$ represents the number of cells disappearing from recruitment in the susceptible compartment and entering into an infected compartment, where f is a nonlinear increasing function.

Using nonlinear incidence rates of the form ${\eta }_{1}f\left(s\right)p$ and ${\eta }_{2}f\left(s\right)y$ into this model is important because the number of effective contacts between infective cells and susceptible cells may increase at high infective levels due to crowding of susceptible cells. If the function $f\left(s\right)$ is increasing for large values of s, used for interpreting the “psychological” effects: for high number of susceptible cells, the infection risk increase as the number of susceptible cells increases, since the total number of cells may tend to increase the number of contacts.

Let $s_0={\displaystyle \frac{\mathrm{\Gamma }}{\delta }}$ and $x_0={\displaystyle \frac{\lambda }{m}}$. We make the following assumption that will be used in the rest of the paper.

Assumption 1.

f is an increasing, nonnegative $C^1\left({\mathbb{R}}_{+}\right)$ function such that $f\left(0\right)=0$ and $f\left(s_0\right)<{\displaystyle \frac{ϵ}{{\eta }_2}}$.

Assumption 1 states that the CHIKV-cell and the infected-cell incidence rates increase with the number of susceptible cells. Assumption 1 also considers that no CHIKV-cell and infected-cell infections can take place in the absence of susceptible cells.

The classical Monod functions can be used to express the transmission rate of infections from infected to susceptible cells as well as from CHIKV to susceptible cells:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle f\left(s\right)={\displaystyle \frac{\overline{\mu }s}{k+s}},}\end{array}\end{array}$$

where $\overline{\mu }$ represents the transmission rates of the disease and k is the Monod constant which is proportional both to the number of CHIKV pathogens and infected cells when the saturated incidence rate is $\overline{\mu }/2$.

The last condition of Assumption 1 is of a mathematical artifice that we used to prove the existence and uniqueness of the infected equilibrium point.

2.1. Basic Properties

First note that $\dot{s}\le \mathrm{\Gamma }-\delta s$; then, one can easily deduce that

Lemma 1.

if $s\left(0\right)\le s_0$, we have $s\left(t\right)\le s_0+\left(s\left(0\right)-s_0\right)e^{-\delta t}\le s_0,\forall t\ge 0$.

The system (1) admits non-negative bounded solutions. In particular,

Lemma 2.

there exist $M_1,M_2,M_3>0$ such that

$$\mathrm{\Omega }=\{(s,y,p,x)\in {\mathbb{R}}_{\ge 0}^4:0\le s,y\le M_1,0\le p\le M_2,0\le x\le M_3\}$$

is a positively invariant bounded set for system (1).

Proof. 

Note that

$$\begin{array}{cc}\hfill \dot{s}{\mid }_{s=0}& =\mathrm{\Gamma }>0,\hfill \\ \hfill \dot{y}{\mid }_{y=0}& ={\eta }_{1}pf\left(s\right)\ge 0,\forall s,p\ge 0,\hfill \\ \hfill \dot{p}{\mid }_{p=0}& =\pi y\ge 0,\forall y\ge 0,\hfill \\ \hfill \dot{x}{\mid }_{x=0}& =\lambda >0.\hfill \end{array}$$

Consider $S_1\left(t\right)=s\left(t\right)+y\left(t\right)$ and $S_2\left(t\right)=p\left(t\right)+{\displaystyle \frac{r}{\rho }}x\left(t\right)$. Then, one has

$$\begin{array}{c}\hfill {\dot{S}}_1\left(t\right)=\mathrm{\Gamma }-\delta s\left(t\right)-ϵy\left(t\right)\le \mathrm{\Gamma }-{\gamma }_1(s\left(t\right)+y\left(t\right))=\mathrm{\Gamma }-{\gamma }_1{S}_1\left(t\right),\end{array}$$

where ${\gamma }_1=min\{\delta ,ϵ\}.$ Hence $S_1\left(t\right)\le M_1,$ if $S_1\left(0\right)\le M_1,$ with $M_1={\displaystyle \frac{\mathrm{\Gamma }}{{\gamma }_1}}$. Therefore, $0\le s\left(t\right),y\left(t\right)\le M_1$ if $0\le s\left(0\right)+y\left(0\right)\le M_1$. Furthermore, one has

$$\begin{array}{cc}\hfill {\dot{S}}_2\left(t\right)& =\pi y\left(t\right)-cp\left(t\right)+{\displaystyle \frac{r}{\rho }}\lambda -{\displaystyle \frac{mr}{\rho }}x\left(t\right)\hfill \\ & \le \pi M_1+{\displaystyle \frac{r}{\rho }}\lambda -{\gamma }_2\left(p\left(t\right)+{\displaystyle \frac{r}{\rho }}x\left(t\right)\right)\hfill \\ & =\pi M_1+{\displaystyle \frac{r}{\rho }}\lambda -{\gamma }_2{S}_2\left(t\right),\hfill \end{array}$$

where, ${\gamma }_2=min\{c,m\}.$ Then $S_2\left(t\right)\le M_2,$ if $S_2\left(0\right)\le M_2,$ with $M_2={\displaystyle \frac{\pi M_1+{\displaystyle \frac{r}{\rho }}\lambda }{{\gamma }_2}}.$ As all variables are nonnegative, therefore, $0\le p\left(t\right)\le M_2$ and $0\le x\left(t\right)\le M_3$ if $0\le p\left(0\right)+{\displaystyle \frac{r}{\rho }}x\left(0\right)\le M_2,$ with $M_3={\displaystyle \frac{\rho M_2}{r}}.$    □

2.2. Basic Reproduction Number and Equilibrium Points

The basic reproduction rate is a dimensionless quantity which allows the measurement of the ability of an infectious cell to spread infection through a given cell’s population immediately after its introduction. From a mathematical point of view, it allows, under certain conditions, the stability of the points of equilibrium of a dynamic system to be established.

Diekmann et al. [27,28] designed a way to calculate the basic reproduction number ${\mathcal{R}}_0$, named the next-generation matrix method, and this was adapted by Van den Driessche and Watmough [29].

Here, $F=\left(\begin{array}{ccc}{\eta }_{2}f\left(s_0\right)& {\eta }_{1}f\left(s_0\right)& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$ and $V=\left(\begin{array}{cccc}ϵ& 0& 0& \\ -\pi & r{x}_0+c& 0\\ 0& -\rho x_0& m\end{array}\right)$. The determinant of V is given by det$\left(V\right)=mϵ(r{x}_0+c)>0$; thus,

$V^{-1}={\displaystyle \frac{1}{mϵ(r{x}_0+c)}}\left(\begin{array}{ccc}m(r{x}_0+c)& 0& 0\\ m\pi & mϵ& 0\\ \pi \rho x_0& ϵ\rho x_0& ϵ(r{x}_0+c)\end{array}\right)$ and the next-generation matrix is given by

$F{V}^{-1}={\displaystyle \frac{1}{mϵ(r{x}_0+c)}}\left(\begin{array}{ccc}m(r{x}_0+c){\eta }_{2}f\left(s_0\right)+m\pi {\eta }_{1}f\left(s_0\right)& mϵ{\eta }_{1}f\left(s_0\right)& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)$.

Then, the basic reproduction number of system (1) is calculated as the spectral radius of the matrix $F{V}^{-1}$:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathcal{R}}_0={\displaystyle \frac{m(r{x}_0+c){\eta }_{2}f\left(s_0\right)+m\pi {\eta }_{1}f\left(s_0\right)}{mϵ(r{x}_0+c)}}={\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}+{\displaystyle \frac{\pi {\eta }_{1}f\left(s_0\right)}{ϵ(r{x}_0+c)}}.}\hfill \end{array}\end{array}$$

Lemma 3.

 

• If ${\mathcal{R}}_0\le 1,$ then (1) admits only $E_0\in \Omega$ as an equilibrium point.
• If ${\mathcal{R}}_0>1,$ then (1) admits two equilibrium points, $E_0\in \Omega$ and $E_1\in \stackrel{\circ }{\Omega }$, and here, $\stackrel{\circ }{\Omega }$ represents the interior of the set Ω.

Proof. 

Let $E(s,y,p,x)$ be an equilibrium point satisfying

$$\begin{array}{cc}\hfill 0& =\mathrm{\Gamma }-\delta s-{\eta }_{1}pf\left(s\right)-{\eta }_{2}yf\left(s\right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill 0& ={\eta }_{1}pf\left(s\right)+{\eta }_{2}yf\left(s\right)-ϵy,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill 0& =\pi y-cp-rxp,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill 0& =\lambda +\rho xp-mx.\hfill \end{array}$$

(5)

The resolution of Equations (2)–(5) gives us the CHIKV-free equilibrium point $E_0=(s_0,0,0,x_0)$.

Moreover, we have

$$\begin{array}{c}\hfill \begin{array}{c}x={\displaystyle \frac{\lambda }{m-\rho p}},\hfill \\ y={\displaystyle \frac{cp+rxp}{\pi }}={\displaystyle \frac{cp}{\pi }}+{\displaystyle \frac{r\lambda p}{\pi (m-\rho p)}},\hfill \\ s={\displaystyle \frac{\mathrm{\Gamma }-ϵy}{\delta }}=s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}},\hfill \\ {\eta }_{1}pf\left(s\right)=(ϵ-{\eta }_{2}f\left(s\right))y.\hfill \end{array}\end{array}$$

We define the function

$$\begin{array}{c}\hfill \begin{array}{ccc}g\left(p\right)\hfill & =\hfill & {\eta }_{1}f\left(s\right)+({\eta }_{2}f\left(s\right)-ϵ){\displaystyle \frac{y}{p}}\hfill \\ & =\hfill & {\eta }_{1}f\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)\hfill \\ & & +\left({\eta }_{2}f\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)-ϵ\right)\left({\displaystyle \frac{c}{\pi }}+{\displaystyle \frac{r\lambda }{\pi (m-\rho p)}}\right).\hfill \end{array}\end{array}$$

Then, we obtain

$$\begin{array}{ccc}\hfill g\left(0\right)& =& {\eta }_{1}f\left(s_0\right)+({\eta }_{2}f\left(s_0\right)-ϵ){\displaystyle \frac{cm+r\lambda }{\pi m}}\hfill \\ & =& ϵ{\displaystyle \frac{cm+r\lambda }{\pi m}}\left({\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}+{\displaystyle \frac{\pi m{\eta }_{1}f\left(s_0\right)}{ϵ(cm+r\lambda )}}-1\right)\hfill \\ & =& ϵ{\displaystyle \frac{cm+r\lambda }{\pi m}}({\mathcal{R}}_0-1)>0if{\mathcal{R}}_0>1.\hfill \end{array}$$

Now, we have $\underset{p\to {(m/\rho )}^{-}}{lim}\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)=-\infty$, and then, $\underset{p\to {(m/\rho )}^{-}}{lim}{\eta }_{1}f\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)<0$ and $\underset{p\to {(m/\rho )}^{-}}{lim}{\eta }_{2}f\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)<0$. One deduces, therefore, that

$$\begin{array}{c}\hfill \underset{p\to {(m/\rho )}^{-}}{lim}g\left(p\right)<0.\end{array}$$

The derivative of the function g is given by

$$\begin{array}{c}\hfill \begin{array}{cc}{\displaystyle g^{\prime }\left(p\right)}\hfill & =-{\eta }_1\left({\displaystyle \frac{cϵ}{\delta \pi }}+{\displaystyle \frac{ϵr\lambda }{\delta \pi }}{\displaystyle \frac{m}{{(m-\rho p)}^2}}\right)f^{\prime }\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)\hfill \\ & -{\eta }_2\left({\displaystyle \frac{c}{\pi }}+{\displaystyle \frac{r\lambda }{\pi (m-\rho p)}}\right)\left({\displaystyle \frac{cϵ}{\delta \pi }}+{\displaystyle \frac{ϵr\lambda }{\delta \pi }}{\displaystyle \frac{m}{{(m-\rho p)}^2}}\right)f^{\prime }\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)\hfill \\ & +{\displaystyle \frac{r\lambda \rho }{\pi {(m-\rho p)}^2}}\left({\eta }_{2}f\left(s_0-{\displaystyle \frac{cϵp}{\pi \delta }}-{\displaystyle \frac{r\lambda ϵp}{\pi \delta (m-\rho p)}}\right)-ϵ\right)\hfill \\ & \le -{\eta }_1\left({\displaystyle \frac{cϵ}{\delta \pi }}+{\displaystyle \frac{ϵr\lambda }{\delta \pi }}{\displaystyle \frac{m}{{(m-\rho p)}^2}}\right)f^{\prime }\left(s\right)\hfill \\ & -{\eta }_2\left({\displaystyle \frac{c}{\pi }}+{\displaystyle \frac{r\lambda }{\pi (m-\rho p)}}\right)\left({\displaystyle \frac{cϵ}{\delta \pi }}+{\displaystyle \frac{ϵr\lambda }{\delta \pi }}{\displaystyle \frac{m}{{(m-\rho p)}^2}}\right)f^{\prime }\left(s\right)+{\displaystyle \frac{r\lambda \rho }{\pi {(m-\rho p)}^2}}\left({\eta }_{2}f\left(s_0\right)-ϵ\right),\hfill \end{array}\end{array}$$

and thus, by Assumption 1,

$$g^{\prime }\left(p\right)\le 0\forall p\in (0,{\displaystyle \frac{m}{\rho }}).$$

Therefore, the equation $g\left(p\right)=0$ admits a unique solution $p_1\in (0,{\displaystyle \frac{m}{\rho }})$. Thus, we obtain

$$\begin{array}{cc}\hfill x_1& ={\displaystyle \frac{\lambda }{m-\rho p_1}},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill y_1& ={\displaystyle \frac{c{p}_1+{\displaystyle \frac{r{p}_1\lambda }{m-\rho p_1}}}{\pi }}={\displaystyle \frac{c{p}_{1}m-\rho c{p}_1^2+r{p}_1\lambda }{\pi (m-\rho p_1)}},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill s_1& =s_0-{\displaystyle \frac{ϵ}{\delta }}{\displaystyle \frac{c{p}_{1}m-\rho c{p}_1^2+r{p}_1\lambda }{\pi (m-\rho p_1)}}\le s_0.\hfill \end{array}$$

(8)

Thus, an infected steady-state $E_1=(s_1,y_1,p_1,x_1)$ exists if ${\mathcal{R}}_0>1$.

Let’s show that $E_1\in \stackrel{\circ }{\mathrm{\Omega }}$. From the steady-state conditions of $E_1$, we have $\mathrm{\Gamma }=\delta s_1+{\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)\to \delta s_1+ϵy_1=\mathrm{\Gamma }\to 0<s_1<{\displaystyle \frac{\mathrm{\Gamma }}{\delta }}\le M_{1}and0<y_1<{\displaystyle \frac{\mathrm{\Gamma }}{ϵ}}\le M_1$.

Using Equations (4) and (5), we obtain

$$\begin{array}{cc}\hfill c{p}_1& =\pi y_1+{\displaystyle \frac{r}{\rho }}\lambda -{\displaystyle \frac{mr}{\rho }}{x}_1\to c{p}_1+{\displaystyle \frac{mr}{\rho }}{x}_1=\pi y_1+{\displaystyle \frac{r}{\rho }}\lambda <\pi M_1+{\displaystyle \frac{r}{\rho }}\lambda ,\hfill \end{array}$$

which means that $p_1<{\displaystyle \frac{\pi M_1+{\displaystyle \frac{r}{\rho }}\lambda }{c}}\le M_{2}and{x}_1<{\displaystyle \frac{\rho }{r}}{\displaystyle \frac{\pi M_1+{\displaystyle \frac{r}{\rho }}\lambda }{m}}\le {\displaystyle \frac{\rho M_2}{r}}=M_3$.

It follows that $E_1\in \stackrel{\circ }{\mathrm{\Omega }}.$    □

2.3. Local Stability

In epidemiology, the basic reproduction number ${\mathcal{R}}_0$ of an infection is the average number of secondary cases caused by an individual with a communicable disease in a fully susceptible population.

More precisely, ${\mathcal{R}}_0$ is the average number of people a contagious person can infect. This rate is calculated from a population that is fully susceptible to infection and which has not yet been vaccinated or immunized against an infectious agent.

Theorem 1.

If ${\mathcal{R}}_0<1$, then the disease-free steady-state $E_0$ is locally asymptotically stable, and if ${\mathcal{R}}_0>1$, it is unstable.

Proof. 

The Jacobian matrix at point $(s,y,p,x)$ is given by:

$$\begin{array}{c}\hfill J=\left(\begin{array}{cccccc}{\displaystyle -\delta -{\eta }_{1}p{f}^{\prime }\left(s\right)-{\eta }_{2}y{f}^{\prime }\left(s\right)}& & -{\eta }_{2}f\left(s\right)& & -{\eta }_{1}f\left(s\right)& 0\\ {\displaystyle {\eta }_{1}p{f}^{\prime }\left(s\right)+{\eta }_{2}y{f}^{\prime }\left(s\right)}& & {\eta }_{2}f\left(s\right)-ϵ& & {\eta }_{1}f\left(s\right)& 0\\ {\displaystyle 0}& & \pi & & -(c+rx)& -rp\\ {\displaystyle 0}& & 0& & \rho x& \rho p-m\end{array}\right).\end{array}$$

Its value at $E_0$ is given by:

$$\begin{array}{c}\hfill J_0=\left(\begin{array}{cccccc}{\displaystyle -\delta }& & -{\eta }_{2}f\left(s_0\right)& & -{\eta }_{1}f\left(s_0\right)& 0\\ {\displaystyle 0}& & {\eta }_{2}f\left(s_0\right)-ϵ& & {\eta }_{1}f\left(s_0\right)& 0\\ {\displaystyle 0}& & \pi & & -(c+r{x}_0)& 0\\ {\displaystyle 0}& & 0& & \rho x_0& -m\end{array}\right).\end{array}$$

$J_0$ admits four eigenvalues; ${\lambda }_1=-\delta <0$ and ${\lambda }_2=-m<0$. ${\lambda }_3$ and ${\lambda }_4$ are eigenvalues of the sub-matrix

$$\begin{array}{c}\hfill S_{j_0}:=\left(\begin{array}{ccc}{\displaystyle {\eta }_{2}f\left(s_0\right)-ϵ}& & {\eta }_{1}f\left(s_0\right)\\ {\displaystyle \pi }& & -(c+r{x}_0)\end{array}\right).\end{array}$$

The trace of $S_{j_0}$ is given by

$$\begin{array}{c}\hfill \begin{array}{ccc}Tr\left(S_{j_0}\right)\hfill & =\hfill & {\eta }_{2}f\left(s_0\right)-ϵ-(c+r{x}_0)\hfill \\ & =\hfill & -(c+r{x}_0)-ϵ\left(1-{\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}\right)\hfill \\ & \le \hfill & -(c+r{x}_0)-ϵ\left(1-{\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}-{\displaystyle \frac{\pi {\eta }_{1}f\left(s_0\right)}{ϵ(c+r{x}_0)}}\right)\hfill \\ & \le \hfill & -(c+r{x}_0)-ϵ\left(1-{\mathcal{R}}_0\right)\hfill \end{array}\end{array}$$

and the determinant of $S_{j_0}$ is given by

$$\begin{array}{c}\hfill \begin{array}{ccc}Det\left(S_{j_0}\right)\hfill & =\hfill & -(c+r{x}_0)\left({\eta }_{2}f\left(s_0\right)-ϵ\right)-\pi {\eta }_{1}f\left(s_0\right)\hfill \\ & =\hfill & -ϵ(c+r{x}_0)\left({\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}-1+{\displaystyle \frac{\pi {\eta }_{1}f\left(s_0\right)}{ϵ(c+r{x}_0)}}\right)\hfill \\ & =\hfill & -ϵ(c+r{x}_0)\left({\mathcal{R}}_0-1\right)\hfill \\ & =\hfill & ϵ(c+r{x}_0)\left(1-{\mathcal{R}}_0\right).\hfill \end{array}\end{array}$$

Then, the steady-state $E_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$, and it is unstable if ${\mathcal{R}}_0>1$.    □

Theorem 2.

If ${\mathcal{R}}_0>1$, then the infected steady-state $E_1$ is locally asymptotically stable.

Proof. 

The Jacobian matrix at a point $E_1=(s_1,y_1,p_1,x_1)$ is given by:

$$\begin{array}{ccc}\hfill J_1& =& \left(\begin{array}{cccccc}{\displaystyle -\delta -{\eta }_1{p}_1{f}^{\prime }\left(s_1\right)-{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)}& & -{\eta }_{2}f\left(s_1\right)& & -{\eta }_{1}f\left(s_1\right)& 0\\ {\displaystyle {\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)}& & {\eta }_{2}f\left(s_1\right)-ϵ& & {\eta }_{1}f\left(s_1\right)& 0\\ {\displaystyle 0}& & \pi & & -(c+r{x}_1)& -r{p}_1\\ 0& & 0& & \rho x_1& \rho p_1-m\end{array}\right).\hfill \end{array}$$

The characteristic polynomial is then given by:

$$\begin{array}{c}\hfill P\left(X\right)=\left|\begin{array}{cccccc}{\displaystyle -X-\delta -{\eta }_1{p}_1{f}^{\prime }\left(s_1\right)-{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)}& & -{\eta }_{2}f\left(s_1\right)& & -{\eta }_{1}f\left(s_1\right)& 0\\ {\displaystyle {\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)}& & -X+{\eta }_{2}f\left(s_1\right)-ϵ& & {\eta }_{1}f\left(s_1\right)& 0\\ {\displaystyle 0}& & \pi & & -X-(c+r{x}_1)& -r{p}_1\\ 0& & 0& & \rho x_1& -X+\rho p_1-m\end{array}\right|\\ \hfill =\left|\begin{array}{cccccc}{\displaystyle -(X+\delta )}& & -(X+ϵ)& & 0& 0\\ {\displaystyle {\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)}& & -X+{\eta }_{2}f\left(s_1\right)-ϵ& & {\eta }_{1}f\left(s_1\right)& 0\\ {\displaystyle 0}& & \pi & & -X-(c+r{x}_1)& -r{p}_1\\ 0& & 0& & \rho x_1& -X+\rho p_1-m\end{array}\right|\\ \hfill =-(X+\delta )\left|\begin{array}{cccc}-X+{\eta }_{2}f\left(s_1\right)-ϵ& & {\eta }_{1}f\left(s_1\right)& 0\\ \pi & & -X-(c+r{x}_1)& -r{p}_1\\ 0& & \rho x_1& -X+\rho p_1-m\end{array}\right|\\ \hfill +(X+ϵ)\left|\begin{array}{cccc}{\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right)& & {\eta }_{1}f\left(s_1\right)& 0\\ 0& & -X-(c+r{x}_1)& -r{p}_1\\ 0& & \rho x_1& -X+\rho p_1-m\end{array}\right|\\ \\ \hfill =(X+\delta )[(X+ϵ-{\eta }_{2}f\left(s_1\right))\left((X+c+r{x}_1)(X+m-\rho p_1)+r\rho p_1{x}_1\right)\\ \hfill -\pi {\eta }_{1}f\left(s_1\right)(X+m-\rho p_1)]+({\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right))(X+ϵ)\\ \hfill \left((X+c+r{x}_1)(X+m-\rho p_1)+r\rho p_1{x}_1\right).\end{array}$$

The characteristic polynomial $P\left(X\right)=0$ if, and only if

$$\begin{array}{c}\hfill \begin{array}{c}[(X+\delta )(X+ϵ-{\eta }_{2}f\left(s_1\right))+({\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right))(X+ϵ)]\\ \left((X+c+r{x}_1)(X+m-\rho p_1)+r\rho p_1{x}_1\right)=\pi {\eta }_{1}f\left(s_1\right)(X+\delta )(X+m-\rho p_1)\end{array}\end{array}$$

or also,

$$\begin{array}{c}\hfill \begin{array}{c}[(X+\delta )(X+ϵ-{\eta }_{2}f\left(s_1\right))+({\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right))(X+ϵ)]=\\ {\displaystyle \frac{\pi {\eta }_{1}f\left(s_1\right)(X+\delta )(X+m-\rho p_1)}{\left(\right(X+c+r{x}_1\left)\right(X+m-\rho p_1)+r\rho p_1{x}_1)}}.\end{array}\end{array}$$

Suppose that X is an eigenvalue with $Re\left(X\right)\ge 0$; then, since $(ϵ-{\eta }_{2}f\left(s_1\right))={\displaystyle \frac{{\eta }_1{p}_{1}f\left(s_1\right)}{y_1}}$ and ${\displaystyle \frac{p_1}{y_1}}={\displaystyle \frac{\pi }{(c+r{x}_1)}}$, the left-hand side satisfies

$$\begin{array}{c}\hfill \begin{array}{ccc}|(X+\delta )(X+ϵ-{\eta }_{2}f\left(s_1\right))+({\eta }_1{p}_1{f}^{\prime }\left(s_1\right)+{\eta }_2{y}_1{f}^{\prime }\left(s_1\right))(X+ϵ)|\hfill & >\hfill & (ϵ-{\eta }_{2}f\left(s_1\right))|X+\delta |\hfill \\ & =\hfill & {\displaystyle \frac{{\eta }_1{p}_{1}f\left(s_1\right)}{y_1}}|X+\delta |\hfill \\ & =\hfill & {\displaystyle \frac{\pi {\eta }_{1}f\left(s_1\right)}{c+r{x}_1}}|X+\delta |,\hfill \end{array}\end{array}$$

and the right-hand side satisfies

$$\begin{array}{c}\hfill \begin{array}{ccc}\left|{\displaystyle \frac{\pi {\eta }_{1}f\left(s_1\right)(X+\delta )(X+m-\rho p_1)}{(X+c+r{x}_1)(X+m-\rho p_1)+r\rho p_1{x}_1}}\right|\hfill & <\hfill & \left|{\displaystyle \frac{\pi {\eta }_{1}f\left(s_1\right)(X+\delta )(X+m-\rho p_1)}{(X+c+r{x}_1)(X+m-\rho p_1)}}\right|\hfill \\ & =\hfill & \pi {\eta }_{1}f\left(s_1\right)\left|{\displaystyle \frac{(X+\delta )}{(X+c+r{x}_1)}}\right|\hfill \\ & \le \hfill & {\displaystyle \frac{\pi {\eta }_{1}f\left(s_1\right)}{c+r{x}_1}}|X+\delta |.\hfill \end{array}\end{array}$$

Hence, the contradiction occurs; thus, $Re\left(X\right)<0$ for all eigenvalues X, and therefore, $E_1$ is locally asymptotically stable.    □

3. Global Stability

Consider the function $G\left(z\right)=z-1-lnz$ that will be used in this section.

Theorem 3.

If ${\mathcal{R}}_0\le 1,$ then the CHIKV-free equilibrium point $E_0$ is globally asymptotically stable.

Proof. 

Assume that ${\mathcal{R}}_0\le 1$ and define the following Lyapunov function $U_0(s,y,p,x)$:

$$U_0(s,y,p,x)=s-s_0-{\int }_{s_0}^s{\displaystyle \frac{f\left(s_0\right)}{f\left(v\right)}}dv+y+{\displaystyle \frac{{\eta }_{1}f\left(s_0\right)}{c+r{x}_0}}\left(p+{\displaystyle \frac{r{x}_0}{\rho }}G\left({\displaystyle \frac{x}{x_0}}\right)\right).$$

Clearly, $U_0(s,y,p,x)>0$ for all $s,y,p,x>0$ and $U_0(s_0,0,0,x_0)=0$. The derivative of $U_0$ with respect to time along the system (1) is given by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{U}}_0=& {\displaystyle \left(1-{\displaystyle \frac{f\left(s_0\right)}{f\left(s\right)}}\right)\left(\Gamma -\delta s-{\eta }_{1}pf\left(s\right)-{\eta }_{2}yf\left(s\right)\right)+{\eta }_{1}pf\left(s\right)+{\eta }_{2}yf\left(s\right)-ϵy}\hfill \\ & {\displaystyle +{\displaystyle \frac{{\eta }_{1}f\left(s_0\right)}{c+r{x}_0}}\left(\pi y-cp-rxp+{\displaystyle \frac{r}{\rho }}(1-{\displaystyle \frac{x_0}{x}})(\lambda +\rho xp-mx)\right)}\hfill \\ \hfill =& {\displaystyle \left(1-{\displaystyle \frac{f\left(s_0\right)}{f\left(s\right)}}\right)(\Gamma -\delta s)+{\eta }_{1}pf\left(s_0\right)+{\eta }_{2}yf\left(s_0\right)-ϵy}\hfill \\ & +{\displaystyle \frac{{\eta }_{1}f\left(s_0\right)}{c+r{x}_0}}\left(\pi y+{\displaystyle \frac{r}{\rho }}(1-{\displaystyle \frac{x_0}{x}})(\lambda -mx)-cp-rxp+r(1-{\displaystyle \frac{x_0}{x}})xp\right)\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \le & {\displaystyle \delta \left(1-{\displaystyle \frac{f\left(s_0\right)}{f\left(s\right)}}\right)(s_0-s)+{\eta }_{1}pf\left(s_0\right)+{\eta }_{2}yf\left(s_0\right)-ϵy}\hfill \\ & +{\displaystyle \frac{{\eta }_{1}f\left(s_0\right)}{c+r{x}_0}}\left(\pi y+{\displaystyle \frac{r}{\rho }}(1-{\displaystyle \frac{x_0}{x}})(\lambda -mx)-p(c+r{x}_0)\right)\hfill \\ \hfill \le & {\displaystyle \delta \left(1-{\displaystyle \frac{f\left(s_0\right)}{f\left(s\right)}}\right)(s_0-s)+{\eta }_{2}yf\left(s_0\right)-ϵy+{\displaystyle \frac{{\eta }_{1}f\left(s_0\right)}{c+r{x}_0}}\left(\pi y+{\displaystyle \frac{r}{\rho }}(1-{\displaystyle \frac{x_0}{x}})(\lambda -mx)\right)}\hfill \\ \hfill \le & {\displaystyle -\delta {\displaystyle \frac{(f\left(s\right)-f\left(s_0\right))}{f\left(s\right)}}(s-s_0)+ϵ\left({\displaystyle \frac{{\eta }_{2}f\left(s_0\right)}{ϵ}}+{\displaystyle \frac{\pi {\eta }_{1}f\left(s_0\right)}{ϵ(c+r{x}_0)}}-1\right)y}\hfill \\ & -{\displaystyle \frac{rm{\eta }_{1}f\left(s_0\right)}{\rho (c+r{x}_0)}}{\displaystyle \frac{{(x-x_0)}^2}{x}}\hfill \\ \hfill \le & {\displaystyle -\delta {\displaystyle \frac{(f\left(s\right)-f\left(s_0\right))}{f\left(s\right)}}(s-s_0)-{\displaystyle \frac{rm{\eta }_{1}f\left(s_0\right)}{\rho (c+r{x}_0)}}{\displaystyle \frac{{(x-x_0)}^2}{x}}+ϵ({\mathcal{R}}_0-1)y.}\hfill \end{array}\end{array}$$

If ${\mathcal{R}}_0\le 1,$ then ${\dot{U}}_0\le 0$ for all $s,y,p,x>0.$ Let $W_0=\{(s,y,p,x):{\dot{U}}_0=0\}$. It can be easily shown that $W_0=\left\{E_0\right\}$. Applying LaSalle’s invariance principle [30] (see [20,21,31,32] for other applications), we deduce that $E_0$ is GAS when ${\mathcal{R}}_0\le 1.$    □

Theorem 4.

For system (1), if ${\mathcal{R}}_0>1,$ then $E_1$ isGASin $\stackrel{\circ }{\mathrm{\Omega }}.$

Proof. 

Let a function $U_1(s,y,p,x)$ be defined as:

$$U_1(s,y,p,x)=s-s_1-{\int }_{s_1}^s{\displaystyle \frac{f\left(s_1\right)}{f\left(v\right)}}dv+y_{1}G\left({\displaystyle \frac{y}{y_1}}\right)+{\displaystyle \frac{{\eta }_1{p}_{1}f\left(s_1\right)}{\pi y_1}}{p}_{1}G\left({\displaystyle \frac{p}{p_1}}\right)+{\displaystyle \frac{r{p}_1{\eta }_{1}f\left(s_1\right)}{\rho \pi y_1}}{x}_{1}G\left({\displaystyle \frac{x}{x_1}}\right).$$

Clearly, $U_1(s,y,p,x)>0$ for all nonnegative variables $s,y,p,x>0$ and $U_1(s_1,y_1,p_1,x_1)$$=0.$ Calculating ${\dot{U}}_1$ along the solutions of the model (1), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{U}}_1=& \left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)\left(\mathrm{\Gamma }-\delta s-{\eta }_{1}pf\left(s\right)-{\eta }_{2}yf\left(s\right)\right)+\left(1-{\displaystyle \frac{y_1}{y}}\right)\left({\eta }_{1}pf\left(s\right)+{\eta }_{2}yf\left(s\right)-ϵy\right)\hfill \\ & +{\displaystyle \frac{{\eta }_1{p}_{1}f\left(s_1\right)}{\pi y_1}}\left(1-{\displaystyle \frac{p_1}{p}}\right)\left(\pi y-cp-rxp\right)+{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}\left(1-{\displaystyle \frac{x_1}{x}}\right)\left(\lambda +\rho xp-mx\right)\hfill \\ \hfill =& \left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)(\mathrm{\Gamma }-\delta s)-{\eta }_{1}pf\left(s\right)-{\eta }_{2}yf\left(s\right)+{\eta }_{1}pf\left(s_1\right)+{\eta }_{2}yf\left(s_1\right)+{\eta }_{1}pf\left(s\right)\hfill \\ & +{\eta }_{2}yf\left(s\right)-ϵy-{\eta }_{1}pf\left(s\right){\displaystyle \frac{y_1}{y}}-{\eta }_2{y}_{1}f\left(s\right)+ϵy_1+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p_{1}y}{p{y}_1}}\hfill \\ & -{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{cp}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{c{p}_1}{\pi y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rxp}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rx{p}_1}{\pi y_1}}\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rxp}{\pi y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_{1}p}{\pi y_1}}+{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}\left(1-{\displaystyle \frac{x_1}{x}}\right)\left(\lambda -mx\right)\hfill \\ \hfill =& \left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)(\mathrm{\Gamma }-\delta s)+{\eta }_{1}pf\left(s_1\right)+{\eta }_{2}yf\left(s_1\right)-ϵy-{\eta }_{1}pf\left(s\right){\displaystyle \frac{y_1}{y}}-{\eta }_2{y}_{1}f\left(s\right)+ϵy_1\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p_{1}y}{p{y}_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{cp}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{c{p}_1}{\pi y_1}}\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rx{p}_1}{\pi y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_{1}p}{\pi y_1}}+{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}\left(1-{\displaystyle \frac{x_1}{x}}\right)\left(\lambda -mx\right)\hfill \end{array}.\end{array}$$

Applying the steady-state conditions for $E_1$: $\mathrm{\Gamma }=\delta s_1+{\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)$, $\epsilon y_1={\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right),c{p}_1+r{x}_1{p}_1=\pi y_1,\lambda +\rho x_1{p}_1=m{x}_1$, we get

$$\begin{array}{c}\hfill \begin{array}{cc}{\dot{U}}_1\hfill & =-\delta {\displaystyle \frac{(s-s_1)(f\left(s\right)-f\left(s_1\right))}{f\left(s\right)}}+\left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)\left({\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)\right)+{\eta }_{1}pf\left(s_1\right)\hfill \\ & +{\eta }_{2}yf\left(s_1\right)-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}-{\eta }_{2}yf\left(s_1\right)-{\eta }_{1}pf\left(s\right){\displaystyle \frac{y_1}{y}}-{\eta }_2{y}_{1}f\left(s\right)+{\eta }_1{p}_{1}f\left(s_1\right)\hfill \\ & +{\eta }_2{y}_{1}f\left(s_1\right)+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p_{1}y}{p{y}_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p(\pi y_1-r{x}_1{p}_1)}{\pi p_1{y}_1}}\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{(\pi y_1-r{x}_1{p}_1)}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rx{p}_1}{\pi y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_{1}p}{\pi y_1}}\hfill \\ & +{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}\left(1-{\displaystyle \frac{x_1}{x}}\right)\left(m{x}_1-\rho x_1{p}_1-mx\right)\hfill \\ & =-\delta {\displaystyle \frac{(s-s_1)(f\left(s\right)-f\left(s_1\right))}{f\left(s\right)}}+\left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)\left({\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)\right)+{\eta }_{1}pf\left(s_1\right)\hfill \\ & -{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}-{\eta }_{1}pf\left(s\right){\displaystyle \frac{y_1}{y}}-{\eta }_2{y}_{1}f\left(s\right)+{\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{y}{y_1}}\hfill \\ & -{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p_{1}y}{p{y}_1}}-{\eta }_{1}pf\left(s_1\right)+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rp{x}_1}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right)-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1{p}_1}{\pi y_1}}\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rx{p}_1}{\pi y_1}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_{1}p}{\pi y_1}}-m{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}{\displaystyle \frac{{(x-x_1)}^2}{x}}-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1{p}_1}{\pi y_1}}\hfill \\ & +{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1^2{p}_1}{\pi x{y}_1}}\hfill \\ & =-\delta {\displaystyle \frac{(s-s_1)(f\left(s\right)-f\left(s_1\right))}{f\left(s\right)}}+\left(1-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}\right)\left({\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)\right)-{\eta }_{1}pf\left(s\right){\displaystyle \frac{y_1}{y}}\hfill \\ & -{\eta }_2{y}_{1}f\left(s\right)+{\eta }_1{p}_{1}f\left(s_1\right)+{\eta }_2{y}_{1}f\left(s_1\right)-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{p_{1}y}{p{y}_1}}+{\eta }_1{p}_{1}f\left(s_1\right)\hfill \\ & -{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1{p}_1}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{rx{p}_1}{\pi y_1}}-m{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}{\displaystyle \frac{{(x-x_1)}^2}{x}}\hfill \\ & -{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1{p}_1}{\pi y_1}}+{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1^2{p}_1}{\pi x{y}_1}}\hfill \\ & =-\delta {\displaystyle \frac{(s-s_1)(f\left(s\right)-f\left(s_1\right))}{f\left(s\right)}}+{\eta }_1{p}_{1}f\left(s_1\right)\left(3-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}-{\displaystyle \frac{p{y}_{1}f\left(s\right)}{p_{1}yf\left(s_1\right)}}-{\displaystyle \frac{p_{1}y}{p{y}_1}}\right)\hfill \\ & +{\eta }_2{y}_{1}f\left(s_1\right)\left(2-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}-{\displaystyle \frac{f\left(s\right)}{f\left(s_1\right)}}\right)-{\eta }_1{p}_{1}f\left(s_1\right){\displaystyle \frac{r{x}_1{p}_1}{\pi y_1}}\left(2-{\displaystyle \frac{x}{x_1}}-{\displaystyle \frac{x_1}{x}}\right)\hfill \\ & -m{\displaystyle \frac{r{\eta }_1{p}_{1}f\left(s_1\right)}{\rho \pi y_1}}{\displaystyle \frac{{(x-x_1)}^2}{x}}\hfill \\ & =-\delta {\displaystyle \frac{(s-s_1)(f\left(s\right)-f\left(s_1\right))}{f\left(s\right)}}-{\displaystyle \frac{{\eta }_{1}f\left(s_1\right)p_1}{\pi y_1}}{\displaystyle \frac{r\lambda }{\rho x_1}}{\displaystyle \frac{{(x-x_1)}^2}{x}}\hfill \\ & +{\eta }_2{y}_{1}f\left(s_1\right)\left(2-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}-{\displaystyle \frac{f\left(s\right)}{f\left(s_1\right)}}\right)+{\eta }_1{p}_{1}f\left(s_1\right)\left(3-{\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}-{\displaystyle \frac{p{y}_{1}f\left(s\right)}{p_{1}yf\left(s_1\right)}}-{\displaystyle \frac{p_{1}y}{p{y}_1}}\right).\hfill \end{array}\end{array}$$

Using the rule

$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{a}_i\ge \sqrt[n]{{\prod }_{i=1}^n{a}_i},$$

(9)

we get ${\displaystyle \frac{1}{2}}\left({\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}+{\displaystyle \frac{f\left(s\right)}{f\left(s_1\right)}}\right)\ge 1$ and ${\displaystyle \frac{1}{3}}\left({\displaystyle \frac{f\left(s_1\right)}{f\left(s\right)}}+{\displaystyle \frac{p{y}_{1}f\left(s\right)}{p_{1}yf\left(s_1\right)}}+{\displaystyle \frac{p_{1}y}{p{y}_1}}\right)\ge 1$. Therefore, ${\dot{U}}_1\le 0$ for all $s,y,p,x>0$ and ${\dot{U}}_1=0$ if, and only if $s=s_1,y=y_1,p=p_1$ and $x=x_1$. Therefore, we deduce that $E_1$ is globally stable by LaSalle’s invariance principle [30] (see [20,21,31,32] for other applications).    □

4. Optimal Susceptible Recruitment Rate

Consider a time-varying rate of recruitment of uninfected cells, $\mathrm{\Gamma }\left(t\right)$ as a control function. Assume further that f is globally Lipschitz with an upper bound ${\displaystyle \overline{f}=\underset{s>0}{sup}f\left(s\right)}$ and a Lipschitz constant L. The control set ${\mathbf{P}}_{\mathbf{ad}}$ is

$${\mathbf{P}}_{\mathbf{ad}}=\{\mathrm{\Gamma }\left(t\right):0\le {\mathrm{\Gamma }}_{min}\le \mathrm{\Gamma }\left(t\right)\le {\mathrm{\Gamma }}_{max},0\le t\le T,\mathrm{\Gamma }\left(t\right)\mathrm{is}\mathrm{Lebesgue}\mathrm{measurable}\}.$$

The aim is to find the optimal values of $\mathrm{\Gamma }\left(t\right)$, $s\left(t\right)$, $y\left(t\right)$, $p\left(t\right)$, and $x\left(t\right)$ that minimize the criterion function:

$$J\left(\mathrm{\Gamma }\right)={\int }_0^T\left(-{\alpha }_{1}s\left(t\right)+{\alpha }_{2}y\left(t\right)+{\displaystyle \frac{{\alpha }_3}{2}}{\mathrm{\Gamma }}^2\left(t\right)\right)dt.$$

For appropriate positive constants ${\alpha }_1,{\alpha }_2$, and ${\alpha }_3$, the aim is to minimize the infected cells and maximize the susceptible ones, while minimizing the control cost.

Since the optimal control problem is linear with respect to the control with bounded states, it is easy to prove the existence of the optimal solution using classical standard results [33].

Existence and Uniqueness of the Solution

Define the variable $\phi ={\left(s,y,p,x\right)}^t$; then, the model (1) takes the simple form

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \dot{\phi }=G\left(\phi \right)=A\phi +F\left(\phi \right),}\hfill \end{array}\end{array}$$

(10)

where $A=\left(\begin{array}{cccc}-\delta & 0& 0& 0\\ 0& -\epsilon & 0& 0\\ 0& \pi & -c& 0\\ 0& 0& 0& -m\end{array}\right)$ and $F\left(\phi \right)=\left(\begin{array}{c}\mathrm{\Gamma }-pf\left(s\right)-yf\left(s\right)\\ pf\left(s\right)+yf\left(s\right)\\ -rxp\\ \lambda +\rho xp\end{array}\right)$.

Proposition 1.

The function G is continuous and uniformly Lipschitz.

Proof. 

The function F is continuous and uniformly Lipschitz, since

$$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \parallel F\left({\phi }_1\right)-F\left({\phi }_2\right){\parallel }_1& =\hfill & {\displaystyle |-p_{1}f\left(s_1\right)-y_{1}f\left(s_1\right)+p_{2}f\left(s_2\right)+y_{2}f\left(s_2\right)|}\hfill \\ & & +|p_{1}f\left(s_1\right)+y_{1}f\left(s_1\right)-p_{2}f\left(s_2\right)-y_{2}f\left(s_2\right)|+r|-x_1{p}_1+x_2{p}_2|+\rho |x_1{p}_1-x_2{p}_2|\hfill \\ & \le \hfill & {\displaystyle 2|p_{1}f\left(s_1\right)+y_{1}f\left(s_1\right)-p_{2}f\left(s_2\right)-y_{2}f\left(s_2\right)|+2max(r,\rho )|x_1{p}_1-x_2{p}_2|}\hfill \\ & \le \hfill & {\displaystyle 2|(p_1+y_1)f\left(s_1\right)-(p_2+y_2)f\left(s_2\right)|+2max(r,\rho )|x_1(p_1-p_2)+p_2(x_1-x_2)|}\hfill \\ & \le \hfill & {\displaystyle 2|(p_1+y_1)f\left(s_1\right)-(p_1+y_1)f\left(s_2\right)+(p_1+y_1)f\left(s_2\right)-(p_2+y_2)f\left(s_2\right)|}\hfill \\ & \hfill & +2max(r,\rho )M_3|p_1-p_2|+2max(r,\rho )M_2|x_1-x_2|\hfill \\ & \le \hfill & {\displaystyle 2|p_1+y_1\left|\right|f\left(s_1\right)-f\left(s_2\right)|+2|f\left(s_2\right)\left|\right|(p_1+y_1)-(p_2+y_2)|}\hfill \\ & \hfill & +2max(r,\rho )M_3|p_1-p_2|+2max(r,\rho )M_2|x_1-x_2|\hfill \\ & \le \hfill & {\displaystyle 2(M_1+M_2)L|s_1-s_2|+2\overline{f}\left(\right|p_1-p_2|+|y_1-y_2\left|\right)}\hfill \\ & \hfill & +2max(r,\rho )M_3|p_1-p_2|+2max(r,\rho )M_2|x_1-x_2|\hfill \\ & \le \hfill & {\displaystyle M\parallel {\phi }_1-{\phi }_2{\parallel }_1,}\hfill \end{array}\end{array}$$

where $M=2max((M_1+M_2)L,\overline{f},\overline{f}+M_{3}max(r,\rho ),M_{2}max(r,\rho ))$. Since

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle \parallel A{\phi }_1-A{\phi }_2{\parallel }_1\le \parallel A{\parallel }_1\parallel {\phi }_1-{\phi }_2{\parallel }_1,}\end{array}\end{array}$$

(11)

where $\mathrm{\parallel }A{\parallel }_1:={sup}_{X\ne 0}{\displaystyle \frac{\parallel AX{\parallel }_1}{\parallel X{\parallel }_1}}$ is the matrix norm of A subordinate to the vector norm $\parallel ·{\parallel }_1$. Therefore,

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle \parallel G\left({\phi }_1\right)-G\left({\phi }_2\right){\parallel }_1\le K\parallel {\phi }_1-{\phi }_2{\parallel }_1,}\end{array}\end{array}$$

(12)

where $K=max(M,\parallel A\parallel )$, and then the function G is uniformly Lipschitz and continuous.    □

Since the function G is uniformly Lipschitz and continuous, then the solution of system (10) exists and is unique.

Let us apply the maximum principle of Pontryagin [33,34,35] to derive the necessary conditions for the considered optimal control and the corresponding states. The Hamiltonian is

$$\begin{array}{cc}\hfill H& =-{\alpha }_{1}s+{\alpha }_{2}y+{\displaystyle \frac{{\alpha }_3}{2}}{\mathrm{\Gamma }}^2+{\lambda }_1\dot{s}+{\lambda }_2\dot{y}+{\lambda }_3\dot{p}+{\lambda }_4\dot{x}\hfill \\ \hfill & =-{\alpha }_{1}s+{\alpha }_{2}y+{\displaystyle \frac{{\alpha }_3}{2}}{\mathrm{\Gamma }}^2+{\lambda }_1(\mathrm{\Gamma }-\delta s-pf\left(s\right)-yf\left(s\right))+{\lambda }_2(pf\left(s\right)+yf\left(s\right)-ϵy)\hfill \\ \hfill & +{\lambda }_3(\pi y-cp-rxp)+{\lambda }_4(\lambda +\rho xp-mx).\hfill \end{array}$$

(13)

For a given optimal control, ${\mathrm{\Gamma }}^{*}$, there exist adjoint functions ${\lambda }_1,{\lambda }_2,{\lambda }_3$ and ${\lambda }_4$ associated to the states $s,y,p$ and x such that:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\dot{\lambda }}_1\hfill & =-{\displaystyle \frac{\partial H}{\partial s}}={\alpha }_1+{\lambda }_1(\delta +{\eta }_{1}p{f}^{\prime }\left(s\right)+{\eta }_{2}y{f}^{\prime }\left(s\right))-{\lambda }_2({\eta }_{1}p{f}^{\prime }\left(s\right)+{\eta }_{2}y{f}^{\prime }\left(s\right)),\hfill \\ {\dot{\lambda }}_2\hfill & =-{\displaystyle \frac{\partial H}{\partial y}}=-{\alpha }_2+{\lambda }_1{\eta }_{2}f\left(s\right)+{\lambda }_2(\epsilon -{\eta }_{2}f\left(s\right))-{\lambda }_3\pi ,\hfill \\ {\dot{\lambda }}_3\hfill & =-{\displaystyle \frac{\partial H}{\partial p}}={\lambda }_1{\eta }_{1}f\left(s\right)-{\lambda }_2{\eta }_{1}f\left(s\right)+{\lambda }_3(c+rx)-{\lambda }_4\rho x,\hfill \\ {\dot{\lambda }}_4\hfill & =-{\displaystyle \frac{\partial H}{\partial x}}={\lambda }_{3}rp+{\lambda }_4(m-\rho p)\hfill \end{array}\right.\end{array}$$

(14)

where ${\lambda }_i\left(T\right)=0$, for $i=1,2,3,4$, are the transversality conditions.

We minimise the Hamiltonian with respect to the control variable at ${\mathrm{\Gamma }}^{*}$. By using the linearity of the Hamiltonian with respect to the control, we check if the optimal control is bang-bang, singular, or a combination. On the interior of the control set, we have

$${\displaystyle \frac{\partial H}{\partial \mathrm{\Gamma }}}={\alpha }_3\mathrm{\Gamma }+{\lambda }_1.$$

(15)

To investigate the singular case, suppose that ${\displaystyle \frac{\partial H}{\partial \mathrm{\Gamma }}}=0$, then the singular value is given by

$${\mathrm{\Gamma }}_{\mathcal{R}msingular}\left(t\right)=-{\displaystyle \frac{{\lambda }_1}{{\alpha }_3}}$$

if

$${\alpha }_3\ne 0\mathrm{and}{\mathrm{\Gamma }}_{min}\le -{\displaystyle \frac{{\lambda }_1}{{\alpha }_3}}\le {\mathrm{\Gamma }}_{max}.$$

By using the bounds for the control $\mathrm{\Gamma }\left(t\right)$, we get:

$$\begin{array}{c}\mathrm{if}{\displaystyle \frac{\partial H}{\partial \mathrm{\Gamma }}}<0\mathrm{at}t,\mathrm{then}{\mathrm{\Gamma }}^{*}\left(t\right)={\mathrm{\Gamma }}_{max},\\ \mathrm{if}{\displaystyle \frac{\partial H}{\partial \mathrm{\Gamma }}}>0\mathrm{at}t,\mathrm{then}{\mathrm{\Gamma }}^{*}\left(t\right)={\mathrm{\Gamma }}_{min},\\ \mathrm{if}{\displaystyle \frac{\partial H}{\partial \mathrm{\Gamma }}}=0,\mathrm{then}{\mathrm{\Gamma }}_{\mathcal{R}msingular}\left(t\right)=-{\displaystyle \frac{{\lambda }_1}{{\alpha }_3}}.\end{array}$$

Hence, the control is optimal at t, provided ${\alpha }_3\ne 0$ and ${\mathrm{\Gamma }}_{min}\le -{\displaystyle \frac{{\lambda }_1}{{\alpha }_3}}\le {\mathrm{\Gamma }}_{max}$.

5. Numerical Results and Conclusions

We used saturated incidence rates of the form ${\eta }_{1}f\left(s\right)p$ and ${\eta }_{2}f\left(s\right)y$ into this model, where f is a nonlinear increasing function. Thus, the incidence rate considered for numerical simulations is nonlinear and of Monod’s type, given by $f\left(s\right)={\displaystyle \frac{s}{k+s}}$, which is globally Lipschitz with Lipschitz constant ${\displaystyle \frac{1}{k}}$ where k is a constant.

5.1. Numerical Results for the Direct Problem

We consider the parameters $k=1,\mathrm{\Gamma }=2,\lambda =0.1,\pi =1,\delta =1,c=1.2,r=1,\rho =0.02$ and $m=0.6$. For $ϵ=0.8,{\eta }_1=1.4$ and ${\eta }_2=2$ then ${\mathcal{R}}_0=4.35>1$, the solution of (1) converges to $E_1$ (Figure 1). This validates the global stability of $E_1=(s_1,y_1,p_1,x_1)$ when ${\mathcal{R}}_0>1$. For $ϵ=1.2,{\eta }_1=0.2$ and ${\eta }_2=0.5$ then ${\mathcal{R}}_0=0.53<1$, the solution of (1) converges to $E_0=(s_0,0,0,x_0)=(2,0,0,0.167)$ (Figure 2).

5.2. Numerical Simulations for the Control Problem

For the numerical resolution of the optimal control problem, we applied an appropriated scheme based on a modified Gauss-Seidel-like finite-difference scheme (see Appendix A). The parameter values were the same as in Figure 1 (where the equilibrium $E_1$ is globally asymptotically stable) but with a variable, $\mathrm{\Gamma }$. Those values were: $\lambda =0.1,\pi =1,\delta =1,c=1.2,ϵ=0.8,r=1,\rho =0.02,m=0.6,{\eta }_1=1.4,{\eta }_2=2$ and $k=1$. Here, $\mathrm{\Gamma }$ is a variable such that $\mathrm{\Gamma }\left(0\right)=0.1$, ${\mathrm{\Gamma }}^{min}=0$ and ${\mathrm{\Gamma }}^{max}=10$.

In Figure 3, Figure 4 and Figure 5, the behaviours of $\mathrm{\Gamma }$ (left), $s\left(t\right),y\left(t\right),p\left(t\right)$ and $x\left(t\right)$ (right) were plotted for several values of ${\alpha }_1$, ${\alpha }_2$, and ${\alpha }_3$. As expected, the number of uninfected cells increased and converged to 0; however, the number of infected cells decreased. Figure 3, Figure 4 and Figure 5 (right) can be compared to Figure 3 in [5], giving the CHIKV pathogenesis.

6. Conclusions

In the present work, we considered and analyzed a mathematical system of differential equations modelling the Chikungunya virus (CHIKV) dynamics by two ways of infection and general incidence rates. The steady-states, which are the CHIKV-free equilibrium point and CHIKV endemic point, were determined and their local and global stability were studied, based on the basic reproduction number (${\mathcal{R}}_0$) obtained by the next-generation matrix. From the analysis, it was found that the disease-free point is locally asymptotically stable if ${\mathcal{R}}_0\le 1$, and the CHIKV endemic point is locally asymptotically stable if ${\mathcal{R}}_0>1$. The global stability of the equilibrium points were carried out using Lyapunov stability, and we confirmed the local stability results. These results are consistent with those obtained when using particular incidence rates as in [18]. We improved and developed the main mathematical techniques used in the previous studies [15,16,17,18,19,20,21]. Furthermore, we obtained the same sharp threshold criteria for the local and global stability of both disease-free and endemic steady-states. Later, we applied an optimal recruitment strategy to minimize the number of infected cells. We thus designed a nonlinear optimal control problem. Some numerical simulations were conducted to visualize the analytical results obtained.

There is still a series of questions of great interest. Current scientific research about the dynamics of a virus infection is mainly focused on solving three major issues: the role of time-delay, the space-time spread of a virus outbreak that occurs in a specific region of the territory, and the role of intrinsic fluctuations.

In this context, we can seek to develop a mathematical model which improves the model presented here by considering three types of infected cells, namely, latently infected cells, short-lived productively infected cells, and long-lived productively infected cells. Such a model can be applied to human immunodeficiency virus (HIV) dynamics.