Modeling Zika Virus Disease Dynamics with Control Strategies

Abstract

In this research, we formulated a fractional-order model for the transmission dynamics of Zika virus, incorporating three control strategies: health education campaigns, the use of insecticides, and preventive measures. We conducted a theoretical analysis of the model, obtaining the disease-free equilibrium and the basic reproduction number, and analyzing the existence and uniqueness of the model. Additionally, we performed model parameter estimation using real data on Zika virus cases reported in Colombia. We found that the fractional-order model provided a better fit to the real data compared to the classical integer-order model. A sensitivity analysis of the basic reproduction number was conducted using computed partial rank correlation coefficients to assess the impact of each parameter on Zika virus transmission. Furthermore, we performed numerical simulations to determine the effect of memory on the spread of Zika virus. The simulation results showed that the order of derivatives significantly impacts the dynamics of the disease. We also assessed the effect of the control strategies through simulations, concluding that the proposed interventions have the potential to significantly reduce the spread of Zika virus in the population.

1. Introduction

Zika virus infection is a vector-borne disease caused by Zika virus, transmitted by mosquitoes of the Flaviviridae family and the genus Flavivirus [1]. The primary mosquito species responsible for transmitting Zika virus are Aedes aegypti and Aedes albopictus, with the virus originally discovered in rhesus monkeys in Uganda’s Zika forest in 1947 [2]. The first human case of Zika virus infection was reported in 1952 [2]. The transmission mechanism of Zika virus is similar to that of other vector-borne diseases such as malaria, trypanosomiasis, dengue, chikungunya, and yellow fever [3]. Specifically, the disease is transmitted to humans when an infected Aedes aegypti mosquito bites a susceptible human during a blood meal. Conversely, a susceptible mosquito acquires the infection when it bites an infected human, although the mosquitoes themselves are not affected by the virus [4]. In humans, the incubation period of Zika virus disease ranges from 3 to 14 days after exposure, and most people infected with the virus (about 1 in 5) do not exhibit symptoms [5]. The most common symptoms of Zika virus disease include fatigue, fever, rash, red eyes, joint and muscle pain, and headaches [6]. Currently, there is no vaccine or specific treatment for Zika virus disease, and management is based on supportive care, such as providing pain relievers [7]. Preventive strategies, including the use of mosquito nets, wearing long-sleeved clothing, applying mosquito repellents, and restricting travel to affected areas, are essential measures for preventing the spread of the disease [8].

Mathematical models have played a crucial role in epidemiology, helping researchers and public health officials understand the spread, dynamics, and control of diseases within populations [9,10,11,12]. In recent years, fractional-order derivatives have gained popularity in disease modeling, as they allow for more flexible and realistic representations of biological systems [10,13,14]. Unlike integer-order derivatives, fractional-order derivatives capture the memory effects and hereditary properties inherent in biological systems [15]. Moreover, it has been observed that the membranes of living organisms’ cells exhibit fractional-order electrical conductance, which can be modeled using fractional-order derivatives [16]. These derivatives often provide a better fit to real epidemiological data [17]. The most commonly used fractional derivatives include Caputo, Riemann–Liouville, and Atangana–Baleanu derivatives [18]. Recent research has developed and analyzed mathematical models incorporating fractional-order derivatives for vector-borne diseases [2,3,19]. For instance, Prasad et al. [2] used the Caputo derivative to study the effect of memory on Zika virus dynamics, demonstrating that fractional-order models offer a more flexible and accurate framework for predicting the virus’s spread. In [3], an optimal control problem using fractional-order derivatives was analyzed, showing a better fit to Zika outbreak data. Similarly, Sharma et al. [19] formulated a model incorporating vaccination as a control strategy, with simulations suggesting that vaccination could eliminate the disease in a population.

Fractional-order models have gained significant attention in recent years for their ability to capture complex dynamics in various epidemiological systems, including malaria and co-infections. One such model is the fractional-order malaria transmission model that incorporates temperature and rainfall dependencies, developed by Gizaw and Deressa [20]. Their work extends traditional models by applying fractional calculus, which provides a more flexible framework for describing the memory and hereditary properties of the disease. Similarly, Attiq ul Rehman et al. [21] proposed a fractional-order malaria model that includes temporary immunity and relapse, further enriching our understanding of malaria dynamics by addressing key biological processes through fractional operators. In another study, Menbiko and Deressa [22] introduced an age-structured malaria model utilizing the Atangana–Baleanu fractional operators, offering a unique perspective on the age-dependent transmission of malaria. Abioye’s research [23] adds to this body of work by exploring a fractional-order model for malaria and COVID-19 co-infection dynamics, revealing the potential interactions between these two infections. Similarly, Kumar [24] extended fractional modeling to study the co-infection dynamics between malaria and filariasis, providing a framework for analyzing simultaneous epidemic outbreaks.

In this study, we formulate a fractional-order model of Zika virus transmission and examine the effects of control strategies. Although various types of fractional derivatives are available, we chose the Caputo derivative for this model. The decision was influenced by several advantages: for a constant function, the Caputo derivative equals zero, consistent with the result for integer-order differential equations. Additionally, the Caputo fractional derivative allows the use of local initial conditions in the model’s formulation [18]. Recent studies have shown that one of the most effective ways to minimize the spread of Zika virus to humans is through prevention measures. Control strategies such as wearing long-sleeved clothing, using physical barriers, and applying insecticide repellents during the day and early evening are crucial to reducing contact between humans and mosquitoes [25]. Health education programs that raise awareness regarding the use of insecticide to eliminate Aedes mosquitoes may also significantly help prevent the spread of Zika virus in the population [7]. To the best of our knowledge, no mathematical model of Zika virus transmission has incorporated the three control strategies mentioned. Therefore, in this study, we propose a fractional-order model that incorporates prevention measures, health education awareness, and the use of insecticides to assess their effects on reducing the spread of Zika virus disease.

The rest of the article is organized as follows: in Section 2, the mathematical model formulation is presented; in Section 3, the basic properties of the model, namely the positivity of variables and the boundedness of trajectories, are presented; the computation of the reproduction number and the existence of model equilibrium are presented in Section 4; the results and discussion are presented in Section 5; and concluding remarks complete the paper in Section 6.

2. Model Formulation

In this study, we employ a fractional-order derivative model based on the Caputo sense to formulate and investigate the transmission dynamics of the Zika virus. The model captures the interaction between Aedes mosquitoes and the human population. The human population is denoted by h, while Aedes mosquitoes are represented by v. The human population is further subdivided into five compartments according to infection status: susceptible ($S_h\left(t\right)$), exposed ($E_h\left(t\right)$), asymptomatic ($A_h\left(t\right)$), infectious ($I_h\left(t\right)$), and recovered ($R_h\left(t\right)$). The total human population at time t is denoted as $N_h\left(t\right)$, and is defined by the following:

$$N_h\left(t\right)=S_h\left(t\right)+E_h\left(t\right)+A_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)$$

(1)

Similarly, the total population of Aedes mosquitoes at time t, denoted by $N_v\left(t\right)$, is divided into three compartments: susceptible ($S_v\left(t\right)$), exposed ($E_v\left(t\right)$), and infectious ($I_v\left(t\right)$). Therefore, the total mosquito population is given by:

$$N_v\left(t\right)=S_v\left(t\right)+E_v\left(t\right)+I_v\left(t\right)$$

(2)

All parameters and variables in the model are assumed to be positive. The parameters are defined as follows: ${\Lambda }_h$ and ${\Lambda }_v$ represent the recruitment rates of humans and Aedes mosquitoes, respectively, assuming all new individuals are susceptible, either through birth or immigration. The natural mortality rates of humans and mosquitoes are denoted by ${\mu }_h$ and ${\mu }_v$, respectively. The rates of progression from the exposed to infectious compartments for humans and mosquitoes are represented by ${\alpha }_h$ and ${\alpha }_v$, respectively. The recovery rates of asymptomatic humans and infectious mosquitoes are given by ${\kappa }_h$ and ${\gamma }_v$, respectively. The transmission rate of infection from infectious mosquitoes to susceptible humans is represented by ${\beta }_h$, while ${\beta }_v$ denotes the transmission rate from infectious humans to susceptible mosquitoes. The reduction in disease transmission from asymptomatic humans to susceptible mosquitoes is represented by ${\tau }_h$, and ${\delta }_v$ denotes the rate at which mosquitoes bite humans. Recent studies suggest that prevention measures are crucial in minimizing the spread of Zika virus [26,27]. Accordingly, the proposed model incorporates three control strategies: the protection of humans from contact with Aedes mosquitoes, represented by ${\eta }_h$; public awareness and education campaigns about the spread and prevention of Zika virus, denoted by ${\theta }_h$; and the use of insecticides to reduce the mosquito population, represented by ${ϵ}_v$. These strategies are essential to preventing the spread of Zika virus in the population. Based on the above assumptions, we have assumed the following flow chart and non-linear differential equations:

From the flowchart diagram in Figure 1, the dynamics of the interaction between humans and mosquitoes is given in patch i by the following system of ordinary differential equations:

$$\left\{{\displaystyle \begin{array}{ccc}{\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{S}_h\left(t\right)}\hfill & =\hfill & {\displaystyle {\Lambda }_h^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\mu }_h^{\varphi }+{\theta }_h^{\varphi })S_h,}\hfill \\ {\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{E}_h\left(t\right)}\hfill & =\hfill & {\displaystyle (1-{\eta }_h^{\varphi }){\delta }_h^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })E_h,}\hfill \\ {\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{A}_h\left(t\right)}\hfill & =\hfill & {\displaystyle (1-{\omega }_h^{\varphi }){\alpha }_h^{\varphi }{E}_h-({\kappa }_h^{\varphi }+{\mu }_h^{\varphi })A_h,}\hfill \\ {\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{I}_h\left(t\right)}\hfill & =\hfill & {\displaystyle {\omega }_h^{\varphi }{\alpha }_h^{\varphi }{E}_h-({\mu }_h^{\varphi }+{\gamma }_h^{\varphi })I_h,}\hfill \\ \stackrel{c}{t_0}{D}_t^{\varphi }{R}_h\left(t\right)\hfill & =\hfill & {\kappa }_h^{\varphi }{A}_h+{\gamma }_h^{\varphi }{I}_h+{\theta }_h^{\varphi }{S}_h-{\mu }_h^{\varphi }{R}_h,\hfill \\ \stackrel{c}{t_0}{D}_t^{\varphi }{S}_v\left(t\right)\hfill & =\hfill & {\Lambda }_v^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_v^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })S_v,\hfill \\ \stackrel{c}{t_0}{D}_t^{\varphi }{E}_v\left(t\right)\hfill & =\hfill & (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_v^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\alpha }_v^{\varphi })E_v,\hfill \\ \stackrel{c}{t_0}{D}_t^{\varphi }{I}_v\left(t\right)\hfill & =\hfill & {\alpha }_h^{\varphi }{E}_h-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })I_v.\hfill \end{array}}\right.$$

(3)

3. Basic Properties

Herein, we study the basic properties of system (3), which are essential in the proof of the stability analysis.

3.1. Positivity of Solutions

In this section, we investigate the asymptotic behavior of orbits starting in the non-negative cone ${\mathbb{R}}_8^{+}$. Obviously, system (3), which is a $C^{\infty }$ differential system, admits a unique maximal solution for any associated Cauchy problem. In what follows, we claim the following theorem:

Theorem 1. 

Let $\left(t_0=0,X_0\left(t_0\right)=\left(S_h\left(0\right),E_h\left(0\right),A_h\left(0\right),I_h\left(0\right),R_h\left(0\right),S_v\left(0\right),E_v\left(0\right),I_v\left(0\right)\right)\right)\in {\mathbb{R}}_{+}^7$ and for $T\in [0,+\infty ]$, $\left([0,T],X\left(t\right)=\left(S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right)\right)\right)$, the maximal solution of the Cauchy problem associated with system (3). Then, for all $t\in [0,T]$, $X\left(t\right)\in {\mathbb{R}}_{+}^8$, $\forall t\in [0,T]$.

Proof. 

Let

$$\Delta =\{\tilde{t}\in [0,T],S_h\left(t\right)>0,E_h\left(t\right)>0,A_h\left(t\right)>0,I_h\left(t\right)>0,R_h\left(t\right)>0,S_v\left(t\right)>0,E_v\left(t\right)>0,I_v\left(t\right)>0\forall t\in [0,\tilde{t}]\}.$$

By the continuity of functions $S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right)$ and $I_v\left(t\right)$, one can see that $\Delta \ne 0$. Let $\tilde{T}=sup\Delta$. Next, we have to show that $\tilde{T}=T$. Suppose $\tilde{T}<T$, then one can suppose that $S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right)$ and $I_v\left(t\right)$ are non-negative on $[0,\tilde{T}]$. At $\tilde{T}$, at least one of the following conditions is satisfied $S_h\left(\tilde{T}\right)=0$, $E_h\left(\tilde{T}\right)=0,A_h\left(\tilde{T}\right)=0,I_h\left(\tilde{T}\right)=0,R_h\left(\tilde{T}\right)=0$, $S_v\left(\tilde{T}\right)=0,E_v\left(\tilde{T}\right)=0$ and $I_v\left(\tilde{T}\right)=0$. Then, from the last equation of system (3), one can obtain the following:

$$\begin{array}{ccc}{\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{I}_v\left(t\right)}\hfill & =\hfill & {\displaystyle {\alpha }_h^{\varphi }{E}_h-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })I_v.}\hfill \end{array}$$

(4)

Integrating Equation (4) from 0 to $\tilde{T}$ yields

$$I_v\left(\tilde{T}\right)=I_v\left(0\right)+{\int }_0^{\tilde{T}}({\alpha }_h^{\varphi }{E}_h-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })I_v)dt>0.$$

Similarly, one can show that $S_h\left(\tilde{T}\right)>0,E_h\left(\tilde{T}\right)>0,A_h\left(\tilde{T}\right)>0,I_h\left(\tilde{T}\right)>0,R_h\left(\tilde{T}\right)>0$, $S_v\left(\tilde{T}\right)>0,$ and $E_v\left(\tilde{T}\right)>0$. This is the contradiction. Then, $\tilde{T}=T$; consequently, the maximal solution $\left(S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right)\right)$ of the Cauchy problem associated with system (3) is positive. This completes the proof. □

3.2. Boundedness of Trajectories

We obtain the following result regarding the boundedness of the trajectories of system (3).

Lemma 1. 

Each non-negative solution of system (3) is bounded on ${\mathbb{R}}_{+}^8$.

Proof. 

First of all, we separate model system (3) into two parts: the human population and the Aedes mosquito population. Adding the first five equations of model system (3) and using Equation (1) yields the following:

$$\begin{array}{c}\hfill {\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }{N}_h\left(t\right)\le {\Lambda }_h-{\mu }_h{N}_h\left(t\right)}\end{array}$$

(5)

By applying the Gronwall Lemma in (5), one can obtain the following:

$$\begin{array}{cc}\hfill N_h\left(t\right)& \le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}}+\left(N_h\left(0\right)-{\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}}\right)exp-{\mu }_h^{\varphi }t.\hfill \end{array}$$

With this in mind, one can suppose that $N_h\left(t\right)\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}}$ for all $t\ge 0$ if $N_{hi}\left(0\right)\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h}}$. This means the following: $S_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h}},E_h\le {\displaystyle \frac{{\Lambda }_h^{\theta }}{{\mu }_h^{\varphi }}},A_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}},I_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{m{u}_h^{\varphi }}}$ and $R_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}}$. Using the same concept, we prove that $N_{vi}\left(t\right)\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_v^{\varphi }}}.$ This means the following: $S_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_v}},E_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_v^{\varphi }}}$ and $,I_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_v^{\varphi }}}$. This completes the proof. □

Corollary 1. 

The region

$$\Omega ={\Omega }_H\times {\Omega }_V,$$

(6)

where

$${\Omega }_H=\left\{(S_h,E_h,A_h,I_h,R_h)\in {\mathbb{R}}_{+}^5,S_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}},E_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}},A_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}},I_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}},R_h\le {\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }}}\right\},$$

(7)

and

$${\Omega }_V=\left\{(S_v,E_v,I_v)\in {\mathbb{R}}_{+}^3,S_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_h^{\varphi }}},E_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_h^{\varphi }}},I_v\le {\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_h^{\varphi }}}\right\},$$

(8)

is invariant and attractive for system (3).

Thus, system (3) is mathematically and ecologically well posed and is sufficient to consider the dynamics of the flow generated by system (3) in $\Omega$.

3.3. Existence and Uniqueness of Solution

In this section, we present the existence and uniqueness of the solution of the model system (3) using the techniques of a fixed-point theorem. First, we denote the Banach space of all continuous real-valued functions equipped with the norm by $\mathcal{B}=\ell \left(\right[0,T],\Re )$, which is defined as follows:

$$\begin{array}{c}\left|\right|S_h,E_h,A_h,I_h,R_h,S_v,E_v,I_v\left|\right|=\left|\right|S_h\left|\right|+\left|\right|E_h\left|\right|+\left|\right|A_h\left|\right|+\left|\right|I_h\left|\right|+\left|\right|R_h\left|\right|+\left|\right|S_v\left|\right|+\left|\right|E_v\left|\right|+\left|\right|I_v\left|\right|,\hfill \end{array}$$

where $\left|\right|S_h\left(t\right)\left|\right|={sup}_{t\in [0,T]}|S_h\left(t\right)|,||E_h\left(t\right)\left|\right|={sup}_{t\in [0,T]}|E_h\left(t\right)|,|A_h\left(t\right)\left|\right|={sup}_{t\in [0,T]}\left|A_h\left(t\right)\right|$, $\left|\right|I_h\left(t\right)\left|\right|={sup}_{t\in [0,T]}|I_h\left(t\right)|,||R_h\left(t\right)\left|\right|={sup}_{t\in [0,T]}|R_h\left(t\right)|,||S_v\left(t\right)\left|\right|={sup}_{t\in [0,T]}\left|S_v\left(t\right)\right|$, $\left|\right|E_v\left(t\right)\left|\right|={sup}_{t\in [0,T]}|E_v\left(t\right)|,||I_v\left(t\right)\left|\right|={sup}_{t\in [0,T]}\left|I_v\left(t\right)\right|$.

We then use the fractional integral operator ${}^{c}I_{0+}^{\varphi }$ on both sides of system (3):

$${\displaystyle {\displaystyle \begin{array}{ccc}{\displaystyle S_h\left(t\right)-S_h\left(0\right)}\hfill & =\hfill & {\displaystyle {}^{c}I_{0+}^{\theta }\left\{{\Lambda }_h^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\mu }_h^{\varphi }+{\theta }_h^{\varphi })S_h\right\},}\hfill \\ {\displaystyle E_h\left(t\right)-E_h\left(0\right)}\hfill & =\hfill & {\displaystyle {}^{c}I_{0+}^{\theta }\left\{(1-{\eta }_h^{\varphi }){\delta }_h^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })E_h\right\},}\hfill \\ A_h\left(t\right)-A_h\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{(1-{\eta }_h^{\varphi }){\delta }_h^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })E_h\right\},\hfill \\ I_h\left(t\right)-I_h\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{{\omega }_h^{\varphi }{\alpha }_h^{\varphi }{E}_h-({\mu }_h^{\varphi }+{\gamma }_h^{\varphi })I_h\right\},\hfill \\ R_h\left(t\right)-R_h\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{{\kappa }_h^{\varphi }{A}_h+{\gamma }_h^{\varphi }{I}_h+{\theta }_h^{\varphi }{S}_h-{\mu }_h^{\varphi }{R}_h\right\},\hfill \\ S_v\left(t\right)-S_v\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{{\Lambda }_v^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })S_v\right\},\hfill \\ E_v\left(t\right)-E_v\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\alpha }_v^{\varphi })E_v\right\},\hfill \\ I_v\left(t\right)-I_v\left(0\right)\hfill & =\hfill & {}^{c}I_{0+}^{\theta }\left\{{\alpha }_h^{\varphi }{E}_h-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })I_v\right\}.\hfill \end{array}}\}}$$

(9)

which implies that, for $k=1,2,3\dots 7$, we have the following:

$$\left.\begin{array}{ccc}{S}_h(t)\hfill & =\hfill & {\displaystyle S_h(0)+\frac{(1-\varphi )}{M(\varphi )}(F_i(t,S_h(t))+\frac{\varphi }{M(\theta )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,S_h(t))d\tau ,}\hfill \\ E_h(t)\hfill & =\hfill & {\displaystyle E_h(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,E_h(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,E_h(t))d\tau ,}\hfill \\ A_h(t)\hfill & =\hfill & {\displaystyle A_h(0)+\frac{(1-\varphi )}{M(\theta )}(F_k(t,E_h(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,E_h(t))d\tau ,}\hfill \\ I_h(t)\hfill & =\hfill & {\displaystyle I_h(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,I_h(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,I_h(t))d\tau ,}\hfill \\ R_h(t)\hfill & =\hfill & {\displaystyle R_h(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,R_h(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,R_h(t))d\tau ,}\hfill \\ S_v(t)\hfill & =\hfill & {\displaystyle S_v(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,S_v(t))+\frac{\theta }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,S_v(t))d\tau ,}\hfill \\ E_v(t)\hfill & =\hfill & {\displaystyle E_v(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,E_v(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,E_v(t))d\tau ,}\hfill \\ I_v(t)\hfill & =\hfill & {\displaystyle I_v(0)+\frac{(1-\varphi )}{M(\varphi )}(F_k(t,I_v(t))+\frac{\varphi }{M(\varphi )}\frac{1}{\Gamma (\varphi )}{\int }_0^t{F}_k(t,I_v(t))d\tau .}\hfill \end{array}\right\}$$

(10)

where

$$\left.\begin{array}{ccc}{F}_k(t,S_h\left(t\right))\hfill & =\hfill & {\Lambda }_h^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\mu }_h^{\varphi }+{\theta }_h^{\varphi })S_h,\hfill \\ F_k(t,E_h\left(t\right))\hfill & =\hfill & (1-{\eta }_h^{\varphi }){\delta }_h^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })E_h,\hfill \\ F_k(t,A_h\left(t\right))\hfill & =\hfill & (1-{\eta }_h^{\varphi }){\delta }_h^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })E_h,\hfill \\ F_k(t,I_h\left(t\right))\hfill & =\hfill & {\displaystyle {\omega }_h^{\varphi }{\alpha }_h^{\varphi }{E}_h-({\mu }_h^{\varphi }+{\gamma }_h^{\varphi })I_h,}\hfill \\ F_k(t,R_h\left(t\right))\hfill & =\hfill & {\displaystyle {\kappa }_h^{\varphi }{A}_h+{\gamma }_h^{\varphi }{I}_h+{\theta }_h^{\varphi }{S}_h-{\mu }_h^{\varphi }{R}_h,}\hfill \\ F_k(t,S_v\left(t\right))\hfill & =\hfill & {\displaystyle {\Lambda }_v^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })S_v,}\hfill \\ F_k(t,E_v\left(t\right))\hfill & =\hfill & {\displaystyle (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }(I_h+{\tau }_h^{\varphi }{A}_h)S_v-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\alpha }_v^{\varphi })E_v,}\hfill \\ F_k(t,I_v\left(t\right))\hfill & =\hfill & {\displaystyle {\alpha }_h^{\varphi }{E}_h-({\mu }_v^{\varphi }+{ϵ}_v^{\varphi })I_v.}\hfill \end{array}\right\}$$

(11)

The kernels N with $0\le Q_k<1$ satisfy the Lipschitz condition in Equation (11) if and only if $S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),$ and $I_v\left(t\right)$ have an upper bound. In general, if supposing that $S_{hi}\left(t\right)$ and $S_{hi}^{\ast }\left(t\right)$ are two functions, we obtain the following:

$$\left.\begin{array}{ccc}\left|\right|F_{k}t,S_h\left(t\right)-F_k(t,S_h^{\ast }\left(t\right))\left|\right|\hfill & =\hfill & \left|\right|{\Lambda }_h^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h{I}_v-({\mu }_h^{\varphi }+{\theta }_h^{\varphi })S_h,\hfill \\ \hfill & \hfill & -\left({\Lambda }_h^{\varphi }-(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h^{\ast }{I}_v-({\mu }_h^{\varphi }+{\theta }_h^{\varphi })S_h^{\ast },\right)\left|\right|,\hfill \\ \hfill & \hfill & =\left({\delta }_v^{\varphi }{\beta }_h^{\varphi }{I}_v+{\mu }_h^{\varphi }+{\theta }_{ij}^{\varphi }\right)\left|\right|S_h-S_h^{\ast }\left|\right|,\hfill \\ \hfill & \hfill & \le \left({\delta }_v^{\varphi }{\beta }_h^{\varphi }\underset{t\in [0,T]}{sup}{I}_h\left(t\right)+{\mu }_h^{\varphi }+{\theta }_h^{\varphi }\right)\left|\right|S_h-S_h^{\ast }\left|\right|,\hfill \\ \hfill & \hfill & =Q_k\left|\right|S_h-S_h^{\ast }\left|\right|.\hfill \end{array}\right\}$$

(12)

where $Q_k={\delta }_v^{\varphi }{\beta }_h^{\varphi }{sup}_{t\in [0,T]}{I}_h\left(t\right)+{\mu }_h^{\varphi }+{\theta }_h^{\varphi }$. Thus:

$$\begin{array}{c}\left|\right|F_{k}t,S_h\left(t\right)-F_k(t,S_h^{\ast }\left(t\right))\left|\right|\le Q_k\left|\right|S_h-S_h^{\ast }\left|\right|\hfill \end{array}$$

(13)

Applying the same techniques as in (12), we obtain the following:

$$\left.\begin{array}{ccc}\left|\right|F_{k}t,E_h\left(t\right)-F_k(t,E_h^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|E_h-E_h^{\ast }\left|\right|,\hfill \\ r\left|\right|F_{k}t,A_h\left(t\right)-F_k(t,A_h^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|A_h-A_h^{\ast }\left|\right|,\hfill \\ r\left|\right|F_{k}t,I_h\left(t\right)-F_k(t,I_h^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|I_h-I_h^{\ast }\left|\right|,\hfill \\ \left|\right|F_{k}t,R_h\left(t\right)-F_k(t,R_h^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|R_h-R_h^{\ast }\left|\right|,\hfill \\ \left|\right|F_{k}t,S_v\left(t\right)-F_k(t,S_v^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|S_v-S_v^{\ast }\left|\right|,\hfill \\ \left|\right|F_{k}t,E_v\left(t\right)-F_k(t,E_v^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|E_v-E_v^{\ast }\left|\right|,\hfill \\ \left|\right|F_{k}t,I_v\left(t\right)-F_k(t,I_v^{\ast }\left(t\right))\left|\right|\hfill & \le \hfill & Q_k\left|\right|I_v-I_v^{\ast }\left|\right|\hfill \end{array}\right\},$$

(14)

whereby $Q_k$ is the Lipschitz constant for the function $F_k(.)$. Indeed, Equation (10) is in recursive form, as follows:

$$\left.\begin{array}{ccc}{\displaystyle S_h\left(t\right)}\hfill & =\hfill & {\displaystyle S_h\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,S_{hi,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,S_{h,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle E_h\left(t\right)}\hfill & =\hfill & {\displaystyle E_h\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,E_{h,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,E_{h,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle A_h\left(t\right)}\hfill & =\hfill & {\displaystyle A_h\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,A_{h,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,E_{h,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle I_h\left(t\right)}\hfill & =\hfill & {\displaystyle I_h\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,I_{h,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,I_{h,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle R_h\left(t\right)}\hfill & =\hfill & {\displaystyle R_h\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,R_{h,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,A_{h,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle S_v\left(t\right)}\hfill & =\hfill & {\displaystyle S_v\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,S_{v,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,S_{v,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle E_v\left(t\right)}\hfill & =\hfill & {\displaystyle E_v\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,E_{v,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,E_{v,n-1}\left(\tau \right))d\tau ,}\hfill \\ {\displaystyle I_v\left(t\right)}\hfill & =\hfill & {\displaystyle I_v\left(0\right)+\frac{(1-\varphi )}{M\left(\varphi \right)}{F}_k(t,E_{vi,n-1}\left(t\right))+\frac{\varphi }{M\left(\varphi \right)}\frac{1}{\Gamma \left(\varphi \right)}{\int }_0^t{(t-\tau )}^{\varphi -1}{F}_k(\tau ,I_{v,n-1}\left(\tau \right))d\tau .}\hfill \end{array}\right\}$$

(15)

Supposing that ${\Psi }_k$, $k=1,2,\dots 8$ is the difference between successive components in (15), we have the following:

$$\left.\begin{array}{ccc}{\Psi }_n^k\hfill & =\hfill & S_{h,n}\left(t\right)-S_{h,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,S_{h,n-1}\left(t\right))-F_k(t,S_{h,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,S_{h,n-1}\left(\tau \right))-F_k(\tau ,S_{h,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & E_{h,n}\left(t\right)-E_{h,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,E_{hi,n-1}\left(t\right))-F_k(t,E_{h,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,E_{h,n-1}\left(\tau \right))-F_k(\tau ,E_{h,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & I_{h,n}\left(t\right)-I_{h,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,I_{h,n-1}\left(t\right))-F_k(t,I_{h,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,I_{h(n-1)}\left(\tau \right))-F_k(\tau ,I_{h,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & R_{h,n}\left(t\right)-R_{h,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,R_{h,n-1}\left(t\right))-F_7(t,R_{h,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,A_{h,n-1}\left(\tau \right))-F_k(\tau ,A_{h,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & S_{v,n}\left(t\right)-S_{v,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,S_{v,n-1}\left(t\right))-F_k(t,S_{v,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,S_{vi,n-1}\left(\tau \right))-F_k(\tau ,S_{vi,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & E_{v,n}\left(t\right)-E_{v,n-1}\left(t\right)={\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}\left(F_k(t,E_{v,n-1}\left(t\right))-F_k(t,E_{v,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,E_{v,n-1}\left(\tau \right))-F_k(\tau ,E_{v,n-2}\left(\tau \right))\right)d\tau ,\hfill \\ {\Psi }_n^k\hfill & =\hfill & I_{v,n}\left(t\right)-I_{v,n-1}\left(t\right)={\displaystyle \frac{1-\theta }{M\left(\theta \right)}}\left(F_k(t,I_{v,n-1}\left(t\right))-F_k(t,I_{v,n-2}\left(t\right))\right)\hfill \\ \hfill & \hfill & +{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{1}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left(F_k(\tau ,I_{v,n-1}\left(\tau \right))-F_k(\tau ,I_{v,n-2}\left(\tau \right))\right)d\tau .\hfill \end{array}\right\}$$

(16)

Considering that $S_{h,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $E_{h,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $A_{h,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $I_{h,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $R_{h,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $S_{v,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right)$, $E_{v,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right),$ and $I_{v,n}\left(t\right)={\displaystyle \sum _{r=1}^n}{\Psi }_r^k\left(t\right).$ Taking the norm on both sides of the equation and using Equations (15) and (16), we have the following:

$$\begin{array}{ccc}\left|\right|{\Psi }_n^k\left(t\right)\left|\right|\hfill & =\hfill & {\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k\left|\right|Q_{n-1}^k\left(t\right)\left|\right|+{\displaystyle \frac{\varphi }{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{\int }_0^t{(t-\tau )}^{\varphi -1}\left|\right|Q_{n-1}^k\left(\tau \right)\left|\right|d\tau ,\hfill \end{array}$$

(17)

In what follows, we state and prove the following theorem based on the results in (17).

Theorem 2. 

The model system (3) has a unique solution for $t\in [0,T]$ if the condition in (19) below is satisfied.

• Let

$$Q_n^k\left(t\right)=\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right),k=1,2,\dots 7,$$

(18)

since the functions $S_h\left(t\right),E_h\left(t\right),A_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),$ and $I_v\left(t\right)$ are bounded and satisfy the Lipschitz condition, and by using (18), we obtain the following:

$$\begin{array}{ccc}\left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|S_{h,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\theta }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|E_{h,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|A_{h,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|I_{h,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|R_{h,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|S_{v,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_{11}}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|E_{v,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n,\hfill \\ \left|\right|Q_n^k\left(t\right)\left|\right|\hfill & \le \hfill & \left|\right|I_{v,n}\left(0\right)\left|\right|{\left({\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}{Q}_k+{\displaystyle \frac{1}{M\left(\varphi \right)}}{\displaystyle \frac{Q_k}{\Gamma \left(\varphi \right)}}{T}^{\varphi }\right)}^n.\hfill \end{array}\}$$

(19)

Therefore, the sequence in (19) exists, and $\left|\right|Q_n^k\left|\right|\to 0$ as $n\to \infty$, $k=1,2,\dots 8$. In addition, when using the triangular inequality in (19), for any s, we have the following:

$$\left.\begin{array}{ccc}\left|\right|S_{h,n+s}-S_{h,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_1^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|E_{h,n+s}-E_{h,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|A_{h,n+s}-A_{h,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|I_{h,n+s}-I_{h,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-k_5}},\hfill \\ \left|\right|R_{h,n+s}-R_{h,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_9^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|S_{v,n+s}-S_{v,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|E_{v,n+s}-E_{v,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}},\hfill \\ \left|\right|I_{v,n+s}-I_{v,n}\left|\right|\hfill & \le \hfill & {\displaystyle \sum _{r=n+1}^{n+s}}{q}_k^r={\displaystyle \frac{q_k^{n+1}-q_k^{n+s+1}}{1-q_k}}.\hfill \end{array}\right\}$$

(20)

where $q_k$, $k=1,2,\dots 8$ are the terms inside the brackets in (19). Thus, $S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),$$R_h\left(t\right),S_v\left(t\right),E_v\left(t\right)$ and $I_v\left(t\right)$ are the Cauchy sequences in $\mathcal{B}$. Therefore, one can note that as $n\to \infty$ in (15), the limit of these sequences is the unique solution of the model system (3). This completes the proofs of the unique solution of the system (3).

4. Basic Reproduction Number and Existence of Equilibria

In this section, we use the next-generation method, as presented in [28], to compute the threshold quantity ${\mathcal{R}}_0$, which determines the persistence and extinction of disease in the population. It is believed that when the basic reproduction number ${\mathcal{R}}_0>1$, the disease persists in the population, and dies out when ${\mathcal{R}}_0<1$. The model system (3) always has a disease-free equilibrium ${\mathcal{E}}^0$ given by

$$\begin{array}{c}\hfill {\mathcal{E}}^0:\left(S_h^0,E_h^0,A_h^0,I_h^0,R_h^0,S_v^0,E_v^0,I_v^0\right)=\left({\displaystyle \frac{{\Lambda }_h^{\varphi }}{{\mu }_h^{\varphi }+{\theta }_h^{\varphi }}},0,0,0,0,{\displaystyle \frac{{\Lambda }_v^{\varphi }}{{\mu }_v^{\varphi }+{ϵ}_v^{\varphi }}},0,0\right).\end{array}$$

Following the next-generation matrix approach method used in [28], the non-negative matrix $\mathcal{F}$ that denotes the generation of new infection and the non-singular matrix $\mathcal{V}$ that denotes the disease transfer among compartments evaluated at ${\mathcal{E}}^0$ are defined as follows:

$${\mathcal{F}}_i=\left(\begin{array}{ccccc}0& 0& 0& 0& (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h^{\ast }\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& (1-{\eta }_h){\delta }_v^{\varphi }{\beta }_v^{\varphi }{\tau }_h^{\varphi }{S}_v^{\ast }& (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_v^{\varphi }{S}_v^{\ast }& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$$

(21)

$${\mathcal{V}}_i=\left(\begin{array}{ccccc}{\alpha }_h^{\varphi }+{\mu }_h^{\varphi }& 0& {\displaystyle 0}& 0& 0\\ -(1-{\omega }_h^{\varphi }){\alpha }_h^{\varphi }& {\mu }_h^{\varphi }+{\kappa }_h^{\varphi }& 0& 0& 0\\ -{\omega }_h^{\varphi }{\alpha }_h^{\varphi }& 0& {\mu }_h^{\varphi }+{\gamma }_h^{\varphi }& 0& 0\\ 0& 0& 0& {\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\kappa }_v^{\varphi }& 0\\ 0& 0& 0& -{\alpha }_v^{\varphi }& {\mu }_v^{\varphi }+{ϵ}_v^{\varphi }\end{array}\right)$$

(22)

Using (21) and (22), the basic reproduction number ${\mathcal{R}}_0$ of the model system (3) is given by the following:

$${\mathcal{R}}_0={\displaystyle \frac{(1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }}{{\mu }_v^{\varphi }+{ϵ}_v^{\varphi }}}\sqrt{{\displaystyle \frac{{\beta }_v^{\varphi }{\beta }_h^{\varphi }{\Lambda }_v^{\varphi }{\Lambda }_h^{\varphi }{\alpha }_v^{\varphi }{\alpha }_h^{\varphi }\left({\omega }_h^{\varphi }({\kappa }_h^{\varphi }+{\mu }_h^{\varphi }-{\tau }_h^{\varphi })+{\tau }_h^{\varphi }\right)}{({\alpha }_h^{\varphi }+{\mu }_h^{\varphi })({\mu }_h^{\varphi }+{\theta }_h^{\varphi })({\kappa }_h^{\varphi }+{\mu }_h^{\varphi })({\mu }_h^{\varphi }+{\gamma }_h^{\varphi })({\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\alpha }_v^{\varphi })}}}$$

The basic reproduction number ${\mathcal{R}}_0$ is defined as the expected number of secondary cases (mosquito or human) produced in a completely susceptible population by one infectious individual (mosquito or human, respectively) during its lifetime. The square root here is due to the fact that the generation of secondary cases in malaria diseases requires two transmission processes. ${\displaystyle \frac{{\alpha }_h^{\varphi }}{{\alpha }_h^{\varphi }+{\mu }_h^{\varphi }}}$ and ${\displaystyle \frac{{\alpha }_v^{\varphi }}{{\alpha }_v^{\varphi }+{\mu }_v^{\varphi }+{ϵ}_v^{\varphi }}}$ are the average life spans of the exposed human and vector populations, respectively.

Theorem 3. 

If ${\mathcal{R}}_0<1,$ then the DFE of system (3) is globally asymptotically stable in Ω; otherwise, it is unstable.

Proof. 

By considering only the infected compartments from (3), one can write the following:

$$\begin{array}{c}\hfill \stackrel{c}{t_0}{D}_t^{\varphi }x=(F-V)x,\end{array}$$

where $x={(E_h,A_h,I_h,E_v,I_v)}^T$, with F and V defined as follows:

$$\mathcal{F}=\left(\begin{array}{ccccc}0& 0& 0& 0& (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_h^{\varphi }{S}_h^{\ast }\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& (1-{\eta }_h){\delta }_v^{\varphi }{\beta }_v^{\varphi }{\tau }_h^{\varphi }{S}_v^{\ast }& (1-{\eta }_h^{\varphi }){\delta }_v^{\varphi }{\beta }_v^{\varphi }{S}_v^{\ast }& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right)$$

(23)

$$\mathcal{V}=\left(\begin{array}{ccccc}{\alpha }_h^{\varphi }+{\mu }_h^{\varphi }& 0& {\displaystyle 0}& 0& 0\\ -(1-{\omega }_h^{\varphi }){\alpha }_h^{\varphi }& {\mu }_h^{\varphi }+{\kappa }_h^{\varphi }& 0& 0& 0\\ -{\omega }_h^{\varphi }{\alpha }_h^{\varphi }& 0& {\mu }_h^{\varphi }+{\gamma }_h^{\varphi }& 0& 0\\ 0& 0& 0& {\mu }_v^{\varphi }+{ϵ}_v^{\varphi }+{\kappa }_v^{\varphi }& 0\\ 0& 0& 0& -{\alpha }_v^{\varphi }& {\mu }_v^{\varphi }+{ϵ}_v^{\varphi }\end{array}\right)$$

(24)

We can observe that by direct calculation, we can verify that $V^{-1}F$ is a non-negative and irreducible matrix and $\rho \left(V^{-1}F\right)={\mathcal{R}}_0.$ It follows from the Perron–Frobenius theorem [29] that $V^{-1}F$ has a positive left eigenvector $\mathbf{w}$ associated with ${\mathcal{R}}_0;$ that is,

$$\begin{array}{c}\hfill \mathbf{w}{V}^{-1}F={\mathcal{R}}_0\mathbf{w}.\end{array}$$

Since $\mathbf{w}{V}^{-1}$ is a positive vector, we propose the following Lyapunov function to study the global stability of disease-free equilibrium:

$$\begin{array}{c}\hfill \mathcal{L}\left(t\right)=\mathbf{w}{V}^{-1}x.\end{array}$$

Differentiating $\mathcal{L}$ using solutions of (3) leads to

$$\begin{array}{ccc}\hfill {\displaystyle \stackrel{c}{t_0}{D}_t^{\varphi }\mathcal{L}\left(t\right)}& =& \mathbf{w}{V}^{-1}\stackrel{c}{t_0}{D}_t^{\varphi }x\le \mathbf{w}{V}^{-1}(F-V)x\hfill \\ & =& ({\mathcal{R}}_0-1)\mathbf{w}x\le 0\mathrm{if}{\mathcal{R}}_0\le 1.\hfill \end{array}$$

It can be easily verified that the largest invariant subset of $\Gamma$ where $\stackrel{c}{t_0}{D}_t^{\varphi }\mathcal{L}\left(t\right)=0$ is the singleton $\left\{{\mathcal{E}}^0\right\}$. Therefore, by LaSalle’s invariance principle [30], ${\mathcal{E}}^0$ is globally asymptotically stable in $\Omega$ when ${\mathcal{R}}_0\le 1$. □

Euler Approximation Scheme of Model (3) Using Caputo Derivative Sense

Here, we discuss an efficient numerical scheme called the Euler fractional approximation method for model (3) with regard to the Atangana–Baleanu–Caputo fractional derivatives used in [13]. The proposed model differential equation equations can be presented in the following form:

$$\left.\begin{array}{c}{\displaystyle {}^{C}D_0^{\varphi }{S}_h\left(t\right)=G_1(t,S_h),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{E}_h\left(t\right)=G_2(t,E_h),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{A}_h\left(t\right)=G_3(t,E_h),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{I}_h\left(t\right)=G_4(t,I_h),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{R}_h\left(t\right)=G_5(t,R_h),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{S}_v\left(t\right)=G_6(t,S_v),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{E}_v\left(t\right)=G_7(t,E_v),}\hfill \\ {\displaystyle {}^{C}D_0^{\varphi }{I}_v\left(t\right)=G_8\left(t{I}_v\right).}\hfill \end{array}\right\}$$

(25)

Next, we present Equation (25) for $G_1(t,S_h)$, which satisfies the Lipschitz condition, where $S_h\left(0\right)$ is the initial conditions. Now, by applying the non-integer operator to Equation (25), one can obtain the following:

$$\begin{array}{c}{\displaystyle S_{hi}\left(t\right)=S_{hi}\left(0\right){+}^C{I}_0^{\alpha }{G}_1(t,S_{hi})}\hfill \end{array}$$

(26)

where ${}^{C}I_0^{\varphi }$ denotes the fractional integer operator with respect to the Caputo derivatives. For the proposed numerical method, we consider an interval length $[0,d]$ with a time step size of $h={\displaystyle \frac{d}{N}}$, where $N\in \mathcal{N}.$ Let $S_{h_k}$ be the numerical approximation of $S_h\left(t\right)$ at $t=t_k,$ where $t_k=0+kh$ and $k=0,123;\dots N.$ Applying the Euler method in (26), we obtain the following Caputo operator formula:

$$\begin{array}{c}{\displaystyle S_{h_{k+1}}\left(t\right)=S_h\left(0\right)+\frac{1-\varphi }{M\left(\varphi \right)}{G}_1\left(S_{h_{k+1}}\right)+\frac{h^{\varphi }}{M\left(\varphi \right)\Gamma \left(\varphi \right)}\sum _{p=0}^k{z}_{k+1,p}{G}_1\left(S_{h{i}_p}\right),k=0,1,2,\dots N-1}\hfill \end{array}$$

(27)

where $z_{k+1,p}={(k+1-p)}^{\varphi }-{(k-p)}^{\varphi },p=0,1,\dots ,k.$ The stability analysis of the proposed model is given by the following theorem:

Theorem 4. 

The numerical approximation scheme (27) is conditionally stable.

Proof. 

Let ${\tilde{S}}_{h_0}$ and ${\tilde{S}}_{h_p}$ be approximations of $S_{h_0}$ and $S_{h_p},p=0,\dots ,k+1.$ From Equation (27), we obtain the following:

$$\begin{array}{c}{S}_{h_{k+1}}+{\tilde{S}}_{h_{k+1}}=S_{h_0}+{\tilde{S}}_{h_0}+{\displaystyle \frac{1-\theta }{M\left(\varphi \right)}}{G}_1(S_{h_{k+1}}+{\tilde{S}}_{h{i}_{k+1}})+{\displaystyle \frac{h^{\varphi }}{M\left(\varphi \right)\Gamma \left(\varphi \right)}}\sum _{p=0}^k{z}_{k+1,p}{G}_1(S_{h_p}+{\tilde{S}}_{h_p})\hfill \end{array}$$

(28)

Using Equation (27) in (28), we obtain the following:

$$\begin{array}{ccc}|{\tilde{S}}_{h_{k+1}}|\hfill & =\hfill & |S_{h_0}+{\displaystyle \frac{1-\varphi }{M\left(\varphi \right)}}[G_1(S_{h_{k+1}}+{\tilde{S}}_{h_{k+1}})-G_1\left(S_{h_{k+1}}\right)]+\hfill \\ \hfill & \hfill & {\displaystyle \frac{\theta h^{\varphi }}{M\left(\varphi \right)\Gamma (\varphi +1)}}\sum _{p=0}^k{z}_{k+1,p}[G_1(S_{h_p}+{\tilde{S}}_{h_p})-G_1\left(S_{h_p}\right)]|\hfill \end{array}$$

(29)

Applying the Lipschitz condition and triangular inequality, one can obtain the following:

$$\begin{array}{c}|{\tilde{S}}_{h_{k+1}}|\le {ϵ}_0+{\displaystyle \frac{(1-\varphi )V_1}{M\left(\varphi \right)}}|{\tilde{S}}_{h_{k+1}}|+{\displaystyle \frac{\varphi h^{\theta }{V}_1}{M\left(\theta \right)\Gamma (\theta +1)}}\sum _{p=0}^k{z}_{k+1,p}|{\tilde{S}}_{h{i}_p}|\hfill \end{array}$$

(30)

where ${ϵ}_0={\max }_{0\le k\le N}\{|{\tilde{S}}_{h{i}_0}|+{\displaystyle \frac{\varphi h^{\varphi }{V}_1{z}_{k,0}}{M\left(\varphi \right)\Gamma (\varphi +1)}}|{\tilde{S}}_{h_0}|\}$. Equation (30) can be further simplified to obtain the following:

$$\begin{array}{c}|{\tilde{S}}_{h_{k+1}}|\le V_1{V}_{1_{\varphi }}{ϵ}_0+{\displaystyle \frac{\varphi h^{\varphi }{V}_1{V}_{1_{\varphi }}}{M\left(\varphi \right)\Gamma (\varphi +1)}}\sum _{p=0}^k{z}_{k+1,p}|{\tilde{S}}_{h{i}_p}|\hfill \end{array}$$

(31)

where $V_{1_{\varphi }}={\displaystyle \frac{M\left(\varphi \right)}{|(\varphi -1)V_1+M\left(\varphi \right)|}}.$ Finally, we obtain $|\widetilde{S_{h_{k+1}}}|\le C{V}_{1_{\varphi }}{ϵ}_0$, and this completes the proof. □

5. Results and Discussion

5.1. Model Parameter Estimations

In this section, we perform model fitting using data about the progression of Zika virus disease from 1 to 36 days in Colombia, as reported in [9], to fit the proposed model (3). Some of the parameters used in the simulations were adopted from the literature, as shown in

Table 1. Definition of model parameters and values.

Symbol	Definition	Value	Units	Source
${\beta }_h$	Disease transmission from mosquito to human	0.00195	${\mathrm{day}}^{-1}$	[14,31]
${\beta }_v$	Disease transmission from human to mosquito	0.63	${\mathrm{day}}^{-1}$	[14,31]
${\mu }_h$	Natural mortality rate of human	$27.9113$	${\mathrm{day}}^{-1}$	[14,31]
${\mu }_v$	Natural mortality rate of vector	$27.9113$	${\mathrm{day}}^{-1}$	[14,31]
${\alpha }_h$	Progression rate of human from incubation to infectious	0.055	${\mathrm{day}}^{-1}$	[3,9]
${\alpha }_v$	Progression rate of vector from incubation to infectious	0.055	${\mathrm{day}}^{-1}$	[3,9]
${\kappa }_h$	Progression rate of human from asymptomatic to recovered class	$0.004$	${\mathrm{day}}^{-1}$	[3,9]
${\gamma }_h$	Progression rate of human from infectious to recovered class	$0.004$	${\mathrm{day}}^{-1}$	[9]
${\Lambda }_h$	New recruitment of human	$287,460$	${\mathrm{day}}^{-1}$	[14,31]
${\Lambda }_v$	New recruitment of Aedes mosquito	$[40,800]$	${\mathrm{day}}^{-1}$	[3]
${\theta }_h$	Human education awareness on Zika virus disease	$[0,1]$	${\mathrm{day}}^{-1}$	fitted
${ϵ}_h$	Rate of use of insecticides	$[0,1]$	${\mathrm{day}}^{-1}$	fitted
${\eta }_h$	Rate of prevention in contact with mosquitoes	$[0,1]$	${\mathrm{day}}^{-1}$	fitted
${\delta }_h$	Rate of mosquito biting on human	$[0,1]$	${\mathrm{day}}^{-1}$	[9]
${\tau }_h$	Reduction factor of disease transmission	$[0,1]$	${\mathrm{day}}^{-1}$	[3,9]
${\omega }_h$	Proportion of human progress to infectious class	$[0,1]$	${\mathrm{day}}^{-1}$	[3]

, and some parameters were estimated using root-mean square error (RMSE) and the following formula:

$$\begin{array}{c}\hfill RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{k=1}^{36}{(I\left(k\right)-\widehat{I}\left(k\right))}^2},\end{array}$$

(32)

where n is the number of weekly reported Zika cases from 1 to 36 weeks. We assumed the initial population as follows: $S_h\left(0\right)=6000,$ $E_h\left(0\right)=50,$ $A_h\left(0\right)=50,$ $I_h\left(0\right)=5,$ $R_h\left(0\right)=10,$ $S_v\left(0\right)=100,$ $E_v\left(0\right)=2,$ and $I_v\left(0\right)=10.$ In addition, from model (3), the generated cases new are obtained using the term $(1-{\eta }_h){\delta }_v{\beta }_v(I_h+{\tau }_h{A}_h)S_v$, which counts the detected cases.

Figure 2 shows the reported cases of Zika virus in Colombia for 36 days, and the simulation result in Figure 3 shows (a) the fitting of the model at $\varphi =0.5$ to the real data of Zika virus infections. From illustrations, we observed that model (3) had a good fit to the reported real cases. In (b), we observe that the fractional-order model had the best fit to the observed data compared to the classical integer model. The numerical simulation in Figure 4 shows the cumulative sum of square error (SSER) for different derivative orders. Overall, one can observe that the minimum estimation error for the reported cases occurs at $\varphi =0.5.$ Figure 5a illustrates the numerical simulation of residuals against the reported real data of Zika virus infection in Colombia. One can observe that the residuals have neither auto-correlation nor partial auto-correlation, suggesting that the model (3) had a good fit to the Zika cases in Colombia reported in [9]. Figure 6 shows the effect of prevention measures on Zika virus transmission. Overall, we noted that as the number of prevention measures increases, there is a reduction in the rate of Zika virus cases reported; hence, the disease decreases in the population.

5.2. Sensitivity Analysis

The numerical simulation results in Figure 7 show the relationship between individual parameters and ${\mathcal{R}}_0$ when all the parameters are simultaneously varied. We utilized the partial rank correlation coefficients, and the results are presented in Figure 7.

The results indicate that increasing specific model parameters can enhance the transmission of the disease. These parameters include the human birth rate (${\Lambda }_h$), which increases the number of susceptible individuals in the population, and the vector birth rate (${\Lambda }_v$), which boosts the number of mosquitoes capable of transmitting the virus. Additionally, increasing the vector biting rate (${\delta }_v$) and the progression rate of vectors from the exposed to infectious class (${\alpha }_v$) further amplifies transmission. Similarly, the disease transmission rate from infected vectors to susceptible humans (${\beta }_h$) and from infected humans to susceptible vectors (${\beta }_v$) also plays a crucial role in accelerating the spread. On the other hand, several parameters can reduce disease transmission when increased. These include the mortality rate of humans (${\mu }_h$) and vectors (${\mu }_v$), as well as control strategies such as health education campaigns (${\theta }_h$), the use of insecticides (${ϵ}_h$), and preventive measures (${\eta }_h$). Preventive measures, which include the use of mosquito nets, wearing long-sleeved clothing, applying mosquito repellents, and restricting travel to affected areas, are essential for preventing the spread of the disease. These findings suggest that increasing public awareness through health education campaigns, along with the use of insecticides and preventive measures, can significantly reduce the potential for disease transmission.

It is worth noting that some of the measures indicated in the PRCC analysis cannot be encouraged as practical strategies for disease control because they may be beyond the influence of policymakers. For example, addressing the rate of human birth or human mortality as a means of controlling Zika virus transmission would not be feasible or ethical.

We further investigated the relationship between ${\mathcal{R}}_0$ and the five model parameters that strongly correlated with it and the results presented in Figure 8. Overall, the results confirm that an increase in prevention measures, including the use of insecticides, health education campaigns, and increasing the vector mortality rate, will reduce the potential for disease transmission; meanwhile, increasing the vector biting rate will increase the potential for disease transmission.

The contour plot in Figure 9a demonstrates that as both ${ϵ}_v$ (use of insecticides) and ${\eta }_h$ (prevention measures) increase, the reproduction number, ${\mathcal{R}}_0$, decreases significantly. This indicates that the combined efforts of health education campaigns and insecticide usage effectively reduce the potential for disease transmission within the community. It is important to note that while these two factors are critical, other variables might also influence the overall dynamics of disease spread. In Figure 9b, we observe a similar trend: as the use of insecticides (${ϵ}_v$) increases and the vector mortality rate (${\delta }_v$) improves, the reproduction number decreases. This suggests that increasing vector mortality, possibly by reducing biting rates through insecticide application, can significantly curb disease transmission. Finally, in Figure 9c, the plot reveals that when preventive measures (${\eta }_h$) are reduced and the vector mortality rate (${\delta }_v$) increases, ${\mathcal{R}}_0$ rises. This implies that despite an increased vector mortality rate, if preventive measures such as health education campaigns are not maintained, the risk of disease transmission within the community can still increase, highlighting the need for an approach that balances vector control and preventive strategies.

5.3. Impact of Memory Effect on the Disease Transmission

To investigate the effect of memory on the spread of Zika virus, we simulated model system (3) for ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_0>1$, as presented in Figure 10 and Figure 11, respectively. The order of derivative $\varphi$ was varied within a reasonable range ($\varphi =0.1,0.3,0.5,0.7$). As mentioned in the literature, when the fractional order $\varphi =1.0$, the fractional-order model based on the Caputo derivative becomes a classical ordinary differential model. Based on this assertion, we observed that as the order of derivatives $\varphi$ reduced from 1, the memory effects of the system increased, and the model solution grew quickly, peaked earlier, and converged to unique equilibrium points. In particular, for ${\mathcal{R}}_0<1$, the model solution converged to the disease-free equilibrium after 10 days, and for ${\mathcal{R}}_0>1$, the model solutions converged to unique endemic equilibrium. Furthermore, the results showed that when the fractional order decreases, the model solution differs. In particular, one can observe that when the fractional order approaches 0, the memory effects become strong and the model solutions converge to their respective equilibrium point earlier than when the order of derivatives approaches 1. This result was also observed in previous studies.

5.4. Effects of Insecticides Use on the Disease Dynamics

To investigate the effects of insecticide use on Zika virus disease, we simulated model system (3) at different values of ${ϵ}_v$, and the other parameter values are presented in Table 1. The results are presented in Figure 12. From the numerical simulation, the results show that the use of insecticides may lead to a reduction in the potential for disease transmission. In particular, one can observe that for ${ϵ}_v=0$, the disease persists in the population; if ${ϵ}_v\ge 0.2$, the disease dies in the population. The results in Figure 12 show that the use of insecticides has a significant effect on public health when aiming to control the spread of disease in the population.

5.5. Effects of Prevention Measures on the Disease Dynamics

To evaluate the effects of prevention measures on disease dynamics, we performed the numerical simulation of the model (3) at different values of ${\eta }_h$. The results are presented in Figure 13. In all scenarios, we observed that the effective use of prevention measures such as wearing long-sleeved clothes and the use of mosquito nets reduce the spread of disease in the population. In particular, one can observe that for ${\eta }_h=0$, the disease persists in the population, and when ${\eta }_h\approx 0$, the model solution converges to the unique disease-free equilibrium point.

5.6. Effects of Health Education Campaigns on the Disease Dynamics

Next, we evaluate the effectiveness of health education campaigns on reducing the spread of Zika virus disease. The effective use of health education campaigns can be implemented by using different media outlets such as newspapers, television, and other social media to report on how people can prevent themselves coming into contact with mosquitoes that transmit Zika virus disease. Here, we simulated model (3) at $\theta =0$ (no effective use of health education campaigns) and $\theta =0.9$ (effective use of health education campaigns), and the results are presented in Figure 14. The results show that for $\theta =0$, the disease persists, and for $\theta =0.9$, the disease dies in the population.

6. Concluding Remarks

This study presents a novel fractional-order mathematical model for Zika virus transmission, which provides valuable insights into the disease dynamics and the potential impact of various public health interventions. By incorporating memory effects through fractional derivatives, we were able to capture the long-term influence of past infections on disease progression and transmission, which is a significant improvement over classical integer-order models. The sensitivity analysis highlighted the critical parameters influencing the basic reproduction number ($R_0$), particularly the importance of preventive measures, vector control, and health education campaigns in reducing disease transmission. Our simulations demonstrate that an increase in these public health strategies can lead to a significant decrease in $R_0$, thus lowering the potential for Zika virus outbreaks. Furthermore, the fractional-order model ($\varphi =0.5$) provided a better fit to real data from the Zika virus outbreak in Colombia, emphasizing the practical relevance of fractional calculus in modeling epidemiological data with memory effects. This shows the importance of fractional-order models in capturing the complexity of vector-borne diseases, where the effects of past states cannot be ignored. These findings have important implications for public health policy, as they suggest that a combination of prevention measures, vector control, and health education can be highly effective in controlling the spread of Zika virus. Furthermore, the model’s ability to incorporate memory effects offers a more accurate and flexible tool for predicting disease dynamics, especially in the context of emerging infectious diseases. While the current model provides substantial insights, future work could explore the inclusion of seasonal variations and migration patterns, both of which can further influence disease dynamics. Incorporating these factors would allow for a more comprehensive analysis of how Zika virus transmission may vary across different environments and populations, providing a stronger foundation for targeted intervention strategies. In summary, this study improves on the application of fractional-order models, demonstrating their efficacy in improving model accuracy and supporting public health decision-making. The results offer a foundation for the further exploration of fractional-order approaches in other vector-borne and infectious diseases.