Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika BIOMATH https://biomath.math.bas.bg/biomath/index.php/biomath B f Biomath Forum ORIGINAL ARTICLE Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Kifah Al-Maqrashi1, Fatma Al-Musalhi2,∗, Ibrahim M. ELmojtaba1, Nasser Al-Salti3 1Department of Mathematics, Sultan Qaboos University, Muscat, Oman kifahhsm@gmail.com elmojtaba@squ.edu.om 0000-0001-6018-4671 2Center for Preparatory Studies, Sultan Qaboos University, Muscat, Oman fatma@squ.edu.om 0000-0003-1108-3879 3Department of Applied Mathematics and Science, National University of Science and Technology, Muscat, Oman alsalti@nu.edu.om 0000-0001-9726-4624 Received: July 8, 2021, Accepted: December 14, 2022, Published: December 22, 2022 Abstract: A mathematical model of Zika virus transmis- sion, incorporating human movement between rural areas and nearby forests, is presented to investigate the role of human movement in the spread of Zika virus infections in human and mosquito populations. Proportions of both susceptible and infected humans living in rural areas are assumed to move to nearby forest areas. Direct, indirect, and vertical transmission routes are incorporated for all populations. A mathematical analysis of the proposed model is presented. The analysis starts with normalizing the proposed model. The positivity and boundedness of solutions to the normalized model are then addressed. The basic reproduction number is calculated using the next- generation matrix method and its relation to the three routes of disease transmission has been presented. The sensitivity analysis of the basic reproduction number to all model parameters is investigated. The analysis also includes the existence and stability of disease-free and endemic equilibrium points. Bifurcation analysis is also carried out. Finally, numerical solutions to the normalized model are obtained to confirm the theoretical results and demonstrate human movement’s role in disease transmis- sion in human and mosquito populations. Keywords: Zika, Vertical Transmission, Basic Repro- duction Number, Stability Analysis, Sensitivity Analysis, Bifurcation Analysis I. INTRODUCTION Zika is an arboviral disease in the genus flavivirus closely related to yellow fever, West Nile (WN), and dengue (DEN) viruses. It was first identified in 1947 in Zika Forest in Uganda during sylvatic yellow fever surveillance in a sentinel rhesus monkey [1]. In 1954, it was reported in humans for the first time in Nigeria [2]. The Zika epidemic was stated as a Public Health Emergency of International Concern (PHEIC) by the World Health Organization (WHO) on February 1st, 2016 [3]. It has attracted global attention since its worldwide spread among tropical and subtropical re- gions. In Yap Island, Micronesia in 2007, the first Zika outbreak occurred among humans [4]. During 2013- 2014 the largest epidemic of Zika ever reported was in French Polynesia [4]. Since 2014, the Zika virus (ZIKV) has continued spreading to other pacific islands [2]. It reached southern and Central America after 2015 and Brazil and the Caribbean were highly affected by ZIKV [4]. Local transmission of ZIKV was realized in 34 countries by March 2016 [5]. ZIKV is transmitted primarily to the human popu- Copyright: © 2022 Kifah Al-Maqrashi, Fatma Al-Musalhi, Ibrahim M. ELmojtaba, Nasser Al-Salti. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. *Corresponding author Citation: Kifah Al-Maqrashi, Fatma Al-Musalhi, Ibrahim M. ELmojtaba, Nasser Al-Salti, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika, Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 1/14 https://biomath.math.bas.bg/biomath/index.php/biomath mailto:kifahhsm@gmail.com mailto:elmojtaba@squ.edu.om https://orcid.org/0000-0001-6018-4671 mailto:fatma@squ.edu.om https://orcid.org/0000-0003-1108-3879 mailto:alsalti@nu.edu.om https://orcid.org/0000-0001-9726-4624 https://creativecommons.org/licenses/by/4.0/ https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika lation by bites of infected female Aedes mosquitoes. Analysts have found 19 species of Aedes mosquitoes competent of carrying Zika infection, but the foremost common is the tropical privateer, Aedes aegypti. The vector (mosquito) can pass into the human population through biting after taking a blood meal from an in- fected human. In addition, sexual interaction, perinatal transmission, and blood transfusion are other routes of spreading ZIKV between humans even months after infection. A pregnant woman can pass Zika to her baby, which can cause genuine birth defects. Infection with Zika increases the chances of the infant developing injury with microcephaly as reported in [6] and Guillian syndrome as reported in [2] from infected mothers [4]. In February 2016, France registered the first sexually transmitted case of ZIKV [3]. Zika disease is characterized by mild symptoms in- cluding fever, headache, maculopapular rash, joint and muscle pain, conjunctivitis, etc. The clinical symptoms duration is within two to seven days after the bites [3]. Most reports show that Zika is a self–limiting febrile disease that could be misidentified as dengue or chikungunya fever [7]. The prevention of mosquito bites and control of vectors by using insecticide, eradication of adult and larval breeding areas is the only possible treatment available till now [8]. Understanding virus transmission and disease epi- demiology through mathematical modelling are of great importance for disease management. Several mathemat- ical models have been developed to study the dynamics and propose control strategies for the transmission of ZIKV disease. In [3], the authors proposed a Zika math- ematical model by assuming the standard incidence type interaction of human-to-human transmission of the illness. Also, they extended their work to include optimal control programs (insecticide-treated bed nets, mosquito-repulsive lotions, and electronic devices) to reduce the biting rate of vectors, and to decline the spread of the disease among the human population. In [8], authors proposed a Zika mathematical model including the applications of prevention, treatments, and insecticide as the best way to minimize the spread of ZIKV disease. In [9], researchers suggested a multi- fold Zika mathematical model. They considered the transmission of the ZIKV in the adult population and infants either directly by vector bites or through verti- cal transmission from mothers. The model shows that asymptomatic individuals magnify the disease weight in the community. It also indicated that postponing conception, coupled with aggressive vector control and personal protection use, decrease the cases of micro- cephaly and transmission of ZIKV. Globally, the survival of around 1.6 billion rustic people depends on products obtained from local forests, in whole or in part. Those individuals live adjacent to the forest and have had simple survival conditions and livelihoods for many generations. They depend on those natural and wild resources to meet their needs [10]. In this paper, a mathematical model of ZIKV is constructed to demonstrate the specific and realistic conditions, where the nearby movement of humans may contribute to the spread of virus infections. This happens when an infected human with mild symptoms, moves from rural areas to nearby forest areas looking for work or food. Additionally, the movement of a susceptible human can affect the spread of infections via contagious mosquitoes in the forest. Hence, in this paper, we have split the vector compartment, based on mosquito location, into rural areas and nearby forest areas. Human movement between rural areas and their interaction with vector populations are illustrated in Figure 1. In addition, sexual and vertical transmissions in the human population are considered. Also, vertical transmission from a contaminated female mosquito to its offspring, is suggested as a component that guaran- tees the upkeep of ZIKV. The paper is organized as follows. The model for- mulation is described in Section 2. The model analysis includes the positivity, boundedness of the solution, basic reproduction number, and sensitivity analysis are discussed in Section 3. Furthermore, stability analysis and bifurcation analysis are presented. A numerical analysis of the model using assumed baseline param- eters is given in Section 4 to illustrate the effects of highly sensitive parameters on the human population. Finally, the conclusion is given in Section 5. II. MODEL DESCRIPTION In this section, we introduce a model for ZIKV trans- mission between humans and vectors in rural areas and nearby forests. We begin the description of the model with the human compartments. We split the human population into susceptible Sh, symptomatic Ih and recovered Rh. Susceptible humans Sh can get infected with Zika via three main routes [11]: via a mosquito bite (vector transmission), via sexual transmission or blood transfusion (direct transmission), or by being passed from a mother to a newborn child (vertical transmission). Zika causes nearly no mortality among humans and has been a public health crisis for a relatively short pe- Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 2/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Fig. 1: Illustrated figure for human movement and their interactions with vector populations of ZIKV. riod of time, so we assume the total human population remains constant: Sh + Ih + Rh = NH. We split the vector population into the rural popula- tion (Sv, Iv) and the nearby forest population (Su, Iu). The overall vector populations at time t are Sv + Iv = NV and Su + Iu = NU . Rural and forest mosquitoes are assumed to be only infected by infectious humans. The infection period of mosquitoes ends when the mosquitoes die. As mosquitoes travel distances of no more than a few kilometres, forest mosquitoes will have a direct interaction only with the human population moving from rural areas to the forest. Hence, we assume that a proportion κ1 of the susceptible individuals may get infected by infectious mosquitoes that live in forests and nearby rural areas Iu due to their movement to forest areas, and a proportion κ2 of the infected individuals are assumed to move from rural areas to the nearby forests such that κ1 > κ2 and hence they may infect mosquitoes that live in forests. The proportion (1 −κ1) of susceptible humans who stay in the rural areas can get infected by infectious mosquitoes that live in rural areas Iv, and a proportion (1 − κ2) of infected individuals who stay in the rural areas may infect mosquitoes that live in rural areas. Moreover, a proportion (1 − κ1) of susceptible in- dividuals can also get the infection by interaction with (1 −κ2) of infectious humans (symptomatic), through sexual transmission or other direct routes. We assume that a fraction ε1 of newborns are affected and enter the symptomatic class. Evidence suggests that the fraction is about 2/3 [12]. We also assume that ZIKV is transmitted vertically in the vector population [13] and this is the main pathway it survives in the colder months. We incorporate vertical transmission ε2,ε3 of the ZIKV in both vector populations, respec- tively. The set of non-linear differential equations that rep- resents the proposed mathematical model is given by: S′h = µHNH −µHε1Ih − (1 −κ1)β1θ1Iv Sh NH −κ1β2θ1Iu Sh NH − (1 −κ1)(1 −κ2)λIh Sh NH −µHSh I′h = µHε1Ih + (1 −κ1)β1θ1Iv Sh NH + κ1β2θ1Iu Sh NH + (1 −κ1)(1 −κ2)λIh Sh NH − (γ + µH)Ih R′h = γIh −µHRh S′v = µV NV −µV ε2Iv − (1 −κ2)β1θ2Sv Ih NH −µV Sv I′v = µV ε2Iv + (1 −κ2)β1θ2Sv Ih NH −µV Iv S′u = µUNU −µUε3Iu −β2θ2κ2Su Ih NH −µUSu I′u = µUε3Iu + β2θ2κ2Su Ih NH −µUIu (1) with non negative initial conditions Sh(0), Ih(0), Rh(0), Sv(0), Iv(0), Su(0), Iu(0). In addition, the parameters and their values of the system are defined in Table I. Let SH = Sh NH , IH = Ih NH , RH = Rh NH , SV = Sv NV , IV = Iv NV , SU = Su NU , IU = Iu NU , such that SH + IH + RH = 1, SV + IV = 1, SU + IU = 1. Thus, the considered model (1) have been normalized and rewritten as follows: S′H = µH −µHε1IH −κ1β1θ1α1IV SH −κ1β2θ1α2IUSH −κ1κ2λIHSH −µHSH I′H = µHε1IH + κ1β1θ1α1IV SH + κ1β2θ1α2IUSH + κ1κ2λIHSH − (γ + µH)IH R′H = γIH −µHRH S′V = µV −µV ε2IV −κ2β1θ2SV IH −µV SV I′V = µV ε2IV + κ2β1θ2SV IH −µV IV S′U = µU −µUε3IU −β2θ2κ2SUIH −µUSU I′U = µUε3IU + β2θ2κ2SUIH −µUIU (2) Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 3/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Fig. 2: Progression diagram of the proposed ZIKV model. Table I: Parameters used in the model (1). Parameter Symbol Value per day Natural death/birth rate of humans µH 1/(68.5*365) [14] Natural death/birth rate of mosquitoes in rural areas µV [0.025-0.125] [15, 16] Natural death/birth rate of mosquito in forest areas µU [0.025-0.125] [15, 16] Biting rate of rural mosquitoes on humans β1 [0.3-1.5] [17] Biting rate of forest mosquitoes on humans β2 [0.3-1.5] [17] Transmission probability from an infectious mosquito to a susceptible human θ1 [0.1–0.75] [17] Transmission probability from an infectious human to a susceptible mosquito θ2 [0.3–0.75] [17] Direct (sexual) transmission rate between humans λ [0.01-0.47] [18] Recovery rate of humans γ [0.07-0.33] [16] Probability of vertical transmission in humans ε1 0.67 [12] Probability of vertical transmission in rural mosquitoes ε2 0.06 [19] Probability of vertical transmission in forest mosquitoes ε3 0.06 [19] Fraction of susceptible humans moving from rural to forest areas κ1 [0-0.5] Fraction of infected humans moving from rural to forest areas κ2 [0-0.5] where κ1 = (1 −κ1), κ2 = (1 −κ2), α1 = NV NH , α2 = NU NH , and with non-negative initial condition X(0) := ( SH(0),IH(0),RH(0), SV (0),IV (0),SU (0),IU (0) )T . III. MODEL ANALYSIS In this section, the positivity of solutions, the positive invariant set, and the basic reproduction number are dis- cussed. Also, sensitivity analysis and results related to stability analysis and bifurcation analysis are presented. A. Positivity of Solutions and Positively Invariant Set It is clear that model (1) together with the given non-negative initial condition has a unique solution. Next, we show that all solutions remain non-negative for all t ∈ [0,∞) for arbitrary choice of initial con- Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 4/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika ditions to have an epidemiological convincing result. The following theorem demonstrates the positivity and boundedness of state variables: Theorem 1. The solutions SH(t), IH(t), RH(t), SV (t), IV (t), SU (t) and IU (t) of system (2) with non-negative initial conditions SH(0), IH(0), RH(0), SV (0), IV (0), SU (0), IU (0) remain non-negative for all time t > 0 in a positively invariant closed set Ω := { (SH,IH,RH,SV ,IV ,SU,IU ) T ∈ R7+ : 0 6 SH(t),IH(t),RH(t), SV (t),IV (t),SU (t),IU (t) 6 1 } . Proof: Assume that the initial conditions of the system (2) are non-negative. Let t1 > 0 be the first time at which there exists at least one component which is equal to zero and other components are non-negative on [0, t1). In the following, we will show that none of the components can be zero at t1. Let’s first assume that SH(t1) = 0 and other components are non-negative on [0, t1). Now, S′H can be written as S′H = µH(1−ε1)+µHε1RH−m1SH−µH(1−ε1)SH, where m1 = κ1β1θ1α1IV + κ1β2θ1α2IU + κ1κ2λIH > 0. Then, at t1, we have SH(t) dt ∣∣∣∣ t=t1 = µH(1 −ε1) + µHε1RH(t1) > 0, which means that SH(t) is strictly monotonically in- creasing at t1, that is SH(t) < SH(t1) for all t ∈ (t1 − ε,t1), where ε > 0. Since SH(t1) = 0, then, SH(t) < 0 on (t1−ε,t1). This leads to a contradiction. Therefore, SH(t) cannot be zero at t1. Now, we assume that IH(t1) = 0 and other compo- nents are non-negative. Then IH(t) dt ∣∣∣∣ t=t1 = SH(t1) ( κ1β1θ1α1IV (t1) + κ1β2θ1α2IU (t1) ) > 0, which means that IH(t) is strictly monotonically in- creasing at t1. Hence, we also get a contradiction. Next, assume that RH(t1) = 0 and other components are non-negative. Then RH(t) dt ∣∣∣∣ t=t1 = γIH(t1) > 0, which again leads to a contradiction. Similarly, one can prove that the remaining com- ponents of vector populations SV (t), IV (t), SU (t), IU (t) cannot be zero at t1. Hence, from the above, we conclude that such a point t1 at which at least one com- ponent is zero does not exist. Hence, all components remain non-negative for all time t > 0. For the positively invariant closed set Ω, we first note that the set Ω is said to be positively invariant if the initial conditions are in Ω implies that( SH(t),IH(t),RH(t),SV (t),IV (t), SU (t),IU (t) )T ∈ Ω. Let Φ(t) = (Φ1(t), Φ2(t), Φ3(t)) T , where Φ1(t) = SH(t) + IH(t) + RH(t), Φ2(t) = SV (t) + IV (t), Φ3(t) = SU (t) + IU (t). Then Φ′(t) =   µH −µHΦ1(t)µV −µV Φ2(t) µU −µU Φ3(t)   . Now, solving for Φ1, Φ2 and Φ3, we get Φ1(t) = 1 − (1 − Φ1(0)) e−µHt, Φ2(t) = 1 − (1 − Φ2(0)) e−µV t, Φ3(t) = 1 − (1 − Φ3(0)) e−µUt, where Φ1(0) = SH(0) + IH(0) + RH(0), Φ2(0) = SV (0) + IV (0), Φ3(0) = SU (0) + IU (0). It is straightforward to conclude that Φ1(t) 6 1 if Φ1(0) 6 1, Φ2(t) 6 1 if Φ2(0) 6 1, Φ3(t) 6 1 if Φ3(0) 6 1. Thus, we have 0 6 SH(t),IH(t),RH(t),SV (t),IV (t), SU (t),IU (t) 6 1 and hence the set Ω is positively invariant set. More- over, the set Ω is a globally attractive set since if Φi(0) > 1 then lim t→∞ Φi(t) = 1 for i = 1, 2, 3. Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 5/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika B. The Basic Reproduction Number The model (2) has a disease-free equilibrium (DFE): Z0 := (S0H, 0, 0,S 0 V , 0,S 0 U, 0) ∈ Ω, where S0H = S 0 V = S 0 U = 1. The number of new infections produced by a typical infected individual in a population at DFE is called the basic reproduction number R0 which can be obtained by applying the Next Generation Method [20]. The next-generation matrix is: PQ−1 =  λκ1κ2 γ + µH(1 −ε1) κ1α1β1θ1 µV (1 −ε2) κ1α2β2θ1 µU (1 −ε3) κ2β1θ2 γ + µH(1 −ε1) 0 0 κ2β2θ2 γ + µH(1 −ε1) 0 0   , where P is the Jacobian of the transmission matrix which describes the production of new infections, whereas Q is the Jacobian of the transition matrix which describes changes in state and they are given by P =  λκ1κ2 κ1α1β1θ1 κ1α2β2θ1κ2β1θ2 0 0 κ2β2θ2 0 0   and Q =  γ + µH(1 −ε1) 0 00 µV (1 −ε2) 0 0 0 µU (1 −ε3)   . The basic reproduction number R0 is the dominant eigenvalue of PQ−1, which can be expressed as: R0 = 1 2 ( RHH + √ R2HH + 4(RHV + RHU ) ) , where RHH = κ1κ2λ γ + µH(1 −ε1) , RHV = α1κ1κ2β 2 1θ1θ2 µV (1 −ε2) ( γ + µH(1 −ε1) ), RHU = κ1κ2α2β 2 2θ1θ2 µU (1 −ε3) ( γ + µH(1 −ε1) ). Note that RHH represents the contribution to the reproduction number due to human-to-human transmis- sion, RHV represents the contribution to the repro- duction number due to interaction between human and vector in a rural area, and RHU represents the con- tribution to the reproduction number due to interaction between human and vector in the forest area. The square root arises from the two “generations” required for an infected vector or host to “reproduce” itself [20]. Moreover, the threshold of the disease occurs at R0 = 1 ⇐⇒ RHH +RHV +RHU = 1. Also, it can be easily proven that R0 < 1 implies RHH + RHV + RHU < 1, which means the disease to die out, all the transmission routes represented by RHH, RHV and RHU need to be reduced. Clearly, this will also imply that R0 > 1 whenever RHH, RHV or RHU is greater than one. C. Sensitivity Analysis of the Basic Reproduction Num- ber A fundamental and valuable numeric value for the study of infectious disease dynamics is the basic re- production number R0, since it predicts whether an outbreak will be expected to continue (when R0 > 1) or die out (when R0 < 1). Sensitivity analysis of the basic reproduction number allows us to determine which model parameters have the most impact on R0. A highly sensitive parameter leads to a high quantitative variation in R0. Moreover, sensitivity analysis highlights the parameters that must be attacked by intervention and treatment strategies. Here, we adopt the elasticity index (normalized forward sensitivity index) [21], ER0P , which computes the relative change of R0 with respect to any parameter P as follows: ER0P = P R0 lim ∆p→0 ∆R0 ∆P = P R0 ∂R0 ∂P . (3) Evidently, the corresponding model parameters will affect the basic reproduction number either positively or negatively. The positive sign of the sensitivity indices of the parameter denotes the increase of the basic reproduction number R0 as that parameter changes, whereas the negative sign of the sensitivity indices of the parameter denotes the decrease of the basic reproduction number R0 as that parameter changes. Moreover, the magnitude denotes the relative impor- tance of the spotlight parameter. The biting rate of mosquitoes is a significant param- eter of the epidemiology of the parasite or pathogen since it is directly related to the basic reproduction rate and affects the dynamics of disease transmission in both areas. It varies depending on the local abundance of vectors, vector host preferences, and host attractiveness [22]. Mosquitoes adjust their preferences depending on the availability of a specific host species to enhance their reproductive success [22]. It is also affected by environmental factors such as temperature, humidity, and larval food sources. The terms (1−κ1)(1−κ2)β21 and κ1κ2β22 in R0 suggest that Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 6/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika the signs of the elasticity indices of the proportion of susceptible humans, moving from rural to forest areas κ1 and the proportion of infected humans moving from rural to forest areas κ2, depend on the biting rates of mosquitoes β1 and β2. To investigate this dependence, we use relation (3) to obtain the following relations which describe the elasticity indices of κ1 and κ2, respectively, in terms of both β1 and β2, taking all other parameters to be fixed, namely, κ1 = 0.3, κ2 = 0.2, α1 = 2, α2 = 3, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 0.07, µU = 0.07, λ = 0.235, γ = 0.16, θ1 = 0.33, θ2 = 0.3: E R0 κ1 = 0.3 ( −0.5874 + 1 4 −1.9325 − 60.1774β22 + 22.5665β21√ 0.6763 + 42.1241β22 + 6.7699β 2 1 ) 0.4112 + 1 2 √ 0.6763 + 42.1241β22 + 6.7699β 2 1 E R0 κ2 = 0.2 ( −0.5140 + 1 4 −1.6909 − 52.6552β22 − 33.8497β21√ 0.6763 + 42.1241β22 + 6.7699β 2 1 ) 0.4112 + 1 2 √ 0.6763 + 42.1241β22 + 6.7699β 2 1 Clearly, the above two equations show that the signs of the sensitivity indices of κ1 and κ2 depend on the values of β1 and β2. This dependence is illustrated in Figure 3 based on the range of values of β1 and β2 given in Table I. Figure 3 shows three different regions for the signs of the indices depending on the values of β1 and β2, namely, a region of negative indices, a region of positive indices, and a region of negative index for κ1 and a positive index for κ2. These three regions are all feasible for low values of the rural mosquito’s biting rate β1 up to a certain limit, taken here to be around β1 = 0.75. In this case, the region of positive indices is only feasible for higher values of the forest mosquito’s biting rate β2, taken here to be higher than 0.8. This means that the movement of susceptible and infected humans from rural to forest areas will have the effect of increasing the basic reproduction number if the disease transmission in the forest areas is relatively higher than in the rural areas. However, for higher values of β1, the region of positive indices is not feasible even for high values of β2. The only feasible positive index is for κ2, whereas the index of κ1 remains negative for all values of β2 since the disease transmission in the rural area is high. Using β1 = 0.35 and β2 = 0.9, which correspond to the region of positive indices, we calculate the elasticity Fig. 3: Signs of the sensitivity indices of the movement rates κ1 and κ2 in terms of the biting rates β1 and β2. Table II: Sensitivity indices and their interpretation. Para. Value Sens. ind. Interpretation β1 0.35 0.36629 β1 by 10% R0 by 37% β2 0.9 0.38926 β2 by 10% R0 by 39% λ 0.235 0.24444 λ by 10% R0 by 24% κ1 0.3 0.01138 κ1 by 10% R0 by 1% κ2 0.2 0.08773 κ2 by 10% R0 by 9% θ1 0.33 0.37778 θ1 by 10% R0 by 38% θ2 0.3 0.37778 θ2 by 10% R0 by 38% γ 0.16 -0.62217 γ by 10% R0 by 62% α1 2 0.18315 α1 by 10% R0 by 18% α2 3 0.19463 α2 by 10% R0 by 19% µH 0.00004 -0.00005 µH by 10% R0 by 0.005% µV 0.07 -0.18315 µV by 10% R0 by 18% µU 0.07 -0.19463 µU by 10% R0 by 19% ε1 0.67 0.000104 ε1 by 10% R0 by 0.01% ε2 0.06 0.01169 ε2 by 10% R0 by 0.11% ε3 0.06 0.01242 ε3 by 10% R0 by 12% indices for all parameters. The obtained values and their interpretations are listed in Table II. Clearly, the most efficacious parameter is the biting rate of forest mosquitoes on humans β2, i.e., it has a strong positive impact on the value of R0. Also, the biting rate of forest mosquitoes on humans β1 has a positive impact on R0. The transmission probabilities per bite – per human θ1 and per mosquito θ2 – have a positive influence on the value of R0. Similarly, one can note that the proportions of movement for susceptible humans κ1 and infected κ2 have a small positive effect on R0, with κ2 having a higher index than κ1. On the Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 7/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika other hand, the recovery rate of humans γ has the most negative sensitivity index, as it will decrease R0 by 62% when it increases by 10%. Clearly, there is a very small positive effect of the vertical transmission of humans and both vectors ε1,ε2,ε3 with ε3 having the highest index among them. Using the parameter values listed in Table II, the ba- sic reproduction number is estimated to be R0 = 2.093. This value is close to the one obtained in [23], in which the authors have estimated from notification data that the basic reproduction number for ZIKV in Rio de Janeiro is R0 = 2.33. D. Local Stability of the DFE Here we discuss the local stability of the DFE by finding the eigenvalues of the linearized system. The following theorem is devoted to the local stability of the DFE, i.e., the disease would be eliminated under certain conditions. Theorem 2. If R0 ≤ 1, the DFE of the model (2) is locally asymptotically stable. If R0 > 1, it is unstable. Proof: The linearized matrix of the system (2) at the disease-free equilibrium Z0 is: JZ0 =  −µH −µHε1 −λκ1κ2 0 0 0 −γ −µH(1 −ε1) + λκ1κ2 0 0 0 γ −µH 0 0 −β1κ2θ2 0 −µV 0 β1κ2θ2 0 0 0 −β2κ2θ2 0 0 0 β2κ2θ2 0 0 −κ1α1β1θ1 0 −κ1α2β2θ1 κ1α1β1θ1 0 κ1α2β2θ1 0 0 0 −µV ε2 0 0 −µV (1 −ε2) 0 0 0 −µU −µUε3 0 0 −µU (1 −ε3)   . It is clear that the system has three negative eigenval- ues which are `1 = −µV , `2 = −µU and `3 = −µH with multiplicity two. The remaining eigenvalues can be found from the characteristic equation k(`) = 0, where k(`) is given by: k(`) = `3 + k1` 2 + k2` + k3, with k1 = ξ(1 −RHH) + µV (1 −ε2) + µU (1 −ε3), k2 = ξµV (1 −ε2)(1 −RHH −RHV ) + ξµU (1 −ε3)(1 −RHH −RHU ) + µV µU (1 −ε2)(1 −ε3), k3 = ξµV µU (1 −ε2)(1 −ε3)(1 −RHH −RHV −RHU ), where ξ = ( γ + µH(1 −ε1) ) . It is clear that k3 > 0 if RHH + RHV + RHU < 1 which also implies that k1 > 0 and k2 > 0. Hence, in order to use Routh’s stability criterion [24] to show that the roots of the above characteristic equation have negative real parts, it remains to show that k1k2 − k3 is positive, that is: k1k2 −k3 = 2ξµV µU (1 −ε2)(1 −ε3)(1 −RHH) + ξ2µV (1 −ε2)(1 −RHH)(1 −RHH −RHV ) + ξ2µU (1 −ε3)(1 −RHH)(1 −RHU −RHV ) + ξµ2V (1 −ε2) 2(1 −RHH −RHV ) + ξµ2U (1 −ε3) 2(1 −RHH −RHU ) + µ2V µU (1 −ε2) 2(1 −ε3) + µV µ 2 U (1 −ε2)(1 −ε3) 2. Clearly, k1k2 −k3 > 0 if and only if RHH +RHV < 1 and RHH + RHU < 1. Therefore, by Routh’s stability criterion, the roots of the characteristic equation k(`) = 0 have negative real parts, and hence we conclude that the DFE is locally asymptotically stable whenever R0 ≤ 1. Otherwise, it is unstable. E. Global Stability of the DFE When the solution of the dynamical system (2) approaches a unique equilibrium point regardless of initial conditions then the equilibrium point is globally asymptotically stable. The global stability of the DFE will ensure that the disease is eliminated under all initial conditions. In this regard, we state and prove the following theorem: Theorem 3. If R0 ≤ 1, the disease-free equilibrium Z0 is globally asymptotically stable on the compact set Ω. Proof: Applying Castillo-Chavez theorem [25], consider the following two compartments: X(t) =   SH(t) RH(t) SV (t) SU (t)   , Y (t) =  IH(t)IV (t) IU (t)   , Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 8/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika which describe the uninfected and infected individuals of the system (2), respectively. So that system (2) can be written as: dX dt = F(X,Y ), dY dt = G(X,Y ), G(X, 0) = 0, where F(X,Y ) and G(X,Y ) are the corresponding right hand side of system (2). To guarantee the global asymptotic stability of the DFE, according to the Castillo-Chavez theorem, the following two conditions must be satisfied: (H1) For dX dt = F(X, 0), X0 = (1, 0, 1, 1)T is globally asymptotically stable. (H2) Ĝ > 0, where Ĝ(X,Y ) = AY −G(X,Y ) and A = DY G(X 0, 0) is an Metzler matrix ∀(X,Y ) ∈ Ω. To check the first condition, we find: F(X, 0) =   −µHSH + µH −µHRH −µV SV + µV −µUSU + µU   . Solving the system of ODEs in (H1), we obtain the following behavior of each component: SH(t) = 1 + SH(0)e −µHt =⇒ lim t→∞ SH(t) = 1, RH(t) = RH(0)e −µHt =⇒ lim t→∞ RH(t) = 0, SV (t) = 1 + SV (0)e −µV t =⇒ lim t→∞ SV (t) = 1, SU (t) = 1 + SU (0)e −µUt =⇒ lim t→∞ SU (t) = 1. Hence, the first condition is satisfied. Now, to check the second condition, we first find: A =  −ξ + λκ1κ2 κ1α1β1θ1 κ1α2β2θ1κ2β1θ2 −µV (1 −ε2) 0 κ2β2θ2 0 −µU (1 −ε3)   , where ξ = ( γ + µH(1 −ε1) ) . Then, Ĝ(X,Y ) = AY −G(X,Y ) = (κ1α1β1θ1IV + κ1α2β2θ1IU + κ1κ2λIH)(1 −SH)β1κ2θ2IH(1 −SV ) β2κ2θ2IH(1 −SU )   Since 0 6 SH 6 1, 0 6 SV 6 1 and 0 6 SU 6 1 then Ĝ > 0 for all (X,Y ) ∈ Ω. Thus, Z0 is globally asymptotically stable provided that R0 ≤ 1. F. Existence of Endemic Equilibrium The endemic equilibrium is the state where the infection cannot be totally eradicated and the disease progression persists in a population at all times but in relatively low frequency. Here, we discuss the existence of endemic equilibrium. Theorem 4. For model (2) there exists an endemic equilibrium Z∗ ∈ Ω whenever R0 > 1. Proof: Let Z∗ := (S∗H,I ∗ H,R ∗ H,S ∗ V ,I ∗ V ,S ∗ U,I ∗ U ) be the endemic equilibrium of the model (2) such that: S∗H = µH − (γ + µH)I∗H µH , R∗H = γI∗H µH , S∗V = µV (1 −ε2) µV (1 −ε2) + κ2β1θ2I∗H , I∗V = κ2β1θ2I ∗ H µV (1 −ε2) + κ2β1θ2I∗H , S∗U = µU (1 −ε3) µU (1 −ε3) + κ2β2θ2I∗H , I∗U = κ2β2θ2I ∗ H µU (1 −ε3) + κ2β2θ2I∗H , and I∗H satisfies the following equation: q1I ∗4 H + q2I ∗3 H + q3I ∗2 H + q4I ∗ H = 0, where q1 = β1β2λθ 2 2κ2κ2(γ + µH), q2 = ξβ1β2κ2κ2θ 2 2µH(1 −RHH) + ξβ2κ2θ2µV (1 −ε2)(γ + µH)(RHH + RHV ) + ξκ2β1θ2µU (1 −ε3)(γ + µH)(RHH + RHU ), q3 = ξβ2κ2θ2µHµV (1 −ε2)(1 −RHH −RHV ) + ξβ1θ2κ2µHµU (1 −ε3)(1 −RHH −RHU ) + ξµV µU (1 −ε2)(1 −ε3)(γ + µH) (RHH + RHV + RHU ), q4 = ξµHµV µU (1 −ε2)(1 −ε3) (1 −RHH −RHV −RHU ), where ξ = ( γ + µH(1 −ε1) ) . Solving the above equation we get I∗H = 0, which corresponds to the DFE (Z0) and the remaining roots satisfy the cubic equation: q1I ∗3 H + q2I ∗2 H + q3I ∗ H + q4 = 0. Clearly, if RHH + RHV + RHU > 1, then the above equation has a positive root, since q1 > 0 and q4 < 0. Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 9/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Fig. 4: Bifurcation figure when λ is taken as a bifurcation parameter of system (2) with a bifurcation value λ∗ = 0.2132 at R0 = 1 and by fixing parameter α1 = 2, α2 = 3, κ1 = 0.3, κ2 = 0.2, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 1/14, µU = 1/14, β1 = 0.15, β2 = 0.1, γ = 0.16, θ1 = 0.33, θ2 = 0.3. Now, note that q3 can be written in terms of q2 as follows: q3 = ξ ( µHθ2 ( β2κ2µV (1 −ε2) + β1κ2µU (1 −ε3) ) + µV µU (1 −ε2)(1 −ε3)(γ + µH) (RHH + RHV + RHU ) ) − µH γ + µH ( q2 − ξβ1β2κ2κ2θ22µH(1 −RHH) ) . To ensure the uniqueness of the positive roots, we apply Descartes’s Sign Rule [26]. There exists a unique positive root when q2 > 0 regardless of the sign of q3 and this happens if RHH < 1 and RHH + RHU + RHV > 1. However, when q2 < 0 and RHH + RHU + RHV > 1 there exist at least one positive root. Note that the existence of three positive roots is only possible when q2 < 0 and q3 > 0. G. Bifurcation Analysis When the stability of a system is changed as a param- eter changes causing the emergence or disappearance of new stable points, then the system is said to undergo bifurcation. In this section, we prove that system (2) has transcritical bifurcation. The proof is based on the Sotomayor theorem described in [27]. Let F be defined as the right-hand side of the system (2) and Z = (SH,IH,RH,SV ,IV ,SU,IU ) T . At R0 = 1, we can check that the constant term of the characteristic equation of JZ0 is zero which implies that JZ0 has a simple zero eigenvalue. Here, we choose λ as a bifurcation parameter such that the bifurcation value corresponding to R0 = 1 is given by: λ∗ = (γ + µH(1 −ε1))(1 −RHV −RHU ) κ1κ2 . Solving J(Z0,λ∗)v = 0, where v = (v1,v2,v3,v4,v5,v6,v7) T is a nonzero right eigenvector of J(Z0,λ∗) corresponding to the zero eigenvalue, we obtain: v =   − γ + µH γ µH γ 1 − β1κ2θ2µH µV γ(1 −ε2) β1κ2θ2µH µV γ(1 −ε2) − β2κ2θ2µH µUγ(1 −ε3) β2κ2θ2µH µUγ(1 −ε3)   v3, v3 6= 0. Next we find the corresponding nonzero left eigenvector w = (w1,w2,w3,w4,w5,w6,w7)T , which satisfies JT (Z0,λ∗) w = 0. We get: w =   0 1 0 0 α1β1κ1θ1 µV (1 −ε2) 0 α2β2κ1θ1 µU (1 −ε3)   w2, w2 6= 0. Model (2) can be written as dZ/dt = F(Z), where F(Z) is the right hand side of the model. Now, we check the conditions of the Sotomayor theorem and begin with finding Fλ(λ∗,Z0): Fλ(λ ∗,Z0) = (0, 0, 0, 0, 0, 0, 0)T . So, the first condition is satisfied: wTFλ(λ ∗,Z0) = 0. Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 10/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Next, we find the Jacobian of Fλ(λ∗,Z) as follows: DFλ(λ ∗,Z) =  −κ1κ2IH −κ1κ2SH 0 0 0 0 0 κ1κ2IH κ1κ2SH 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0   . Checking the second condition, we have: wTDFλ(λ ∗,Z0)v = κ1κ2w2v2 6= 0. Finally, we check the third condition by finding D2F(λ∗,Z0), where D2 denotes the matrix of the partial derivatives of each component of DF(Z) and we get: D2F(λ∗,Z0)(v,v) =  −2λκ1κ2v1v2 − 2α1κ1β1θ1v1v5 − 2κ1α2β2θ1v1v7 2λκ1κ2v1v2 + 2α1κ1β1θ1v1v5 + 2κ1α2β2θ1v1v7 0 −2κ2β1θ2v4v2 2κ2β1θ2v4v2 −2κ2β2θ2v6v2 2κ2β2θ2v6v2   Thus, wT ( D2F(λ∗,Z0)(v,v) ) = 2v1(α2β2θ1κ1v7 + α1β1θ1κ1v5 + λκ1κ2) + 2β1θ2κ2v2v4w5 + 2β2θ2κ2v2v6w7. By substituting the values of v’s and w’s, we get: wT ( D2F(λ∗,Z0)(v,v) ) = − ( 2µH(γ + µH)(γ + µH(1 −ε1)) γ2 + 2β32κ1κ 2 2θ 2 2θ1α2 γ2µ2U (1 −ε3)2 ) w2v 2 3, which is nonzero since w2 and v3 are nonzero. Hence, the system (2) experiences a transcritical bifurcation at Z0 as the parameter λ passes through the bifurcation value λ = λ∗. The bifurcation diagram is created by using the MATCONT package [28] and is illustrated in Figure 4. This leads us to establish the following theorem: Theorem 5. Model (2) undergoes transcritical bifur- cation at the DFE (Z0) when the parameter λ passes through the bifurcation value λ = λ∗. Remark. We can establish the local stability of en- demic equilibrium using the above calculations. We note that based on Theorem 4 in [20], a and b are given by: a = 1 2 wT ( D2Fλ(λ ∗,Z0)(v,v) ) = 1 2 n∑ i,j,k=1 vivjwk ∂2Fi ∂xj∂xk (λ∗,Z0), b = wTDFλ(λ ∗,Z0)v = n∑ i,j=1 viwj ∂2Fi ∂xj∂λ (λ∗,Z0). According to the calculations in this section, it is clear that b 6= 0 and a < 0 if w2 is positive. Thus, there exists δ > 0 such that the endemic equilibrium Z∗ is locally asymptotically stable near Z0 for 0 < λ < δ. Moreover, according to Castillo-Chavez and Song [29] the direction of the bifurcation of the system (2) at R0 = 1 is forward (supercritical bifurcation). IV. NUMERICAL ANALYSIS In this section, the forgoing theoretical results are confirmed by presenting the numerical results of the Zika SIR-SI model (2). The asymptotic behavior of the model is characterized by solving the system numeri- cally using the numerical simulations of MATLAB with baseline parameters listed in Table I with appropriate initial conditions. We assume that the human population size is 200000, the rural mosquito population size is 400000, and the forest mosquito population size is 600000. These values are obtained by taking into consideration desired conditions or from literature. Phase diagram for the case R0 < 1 is illustrated in Figure 5. Here, the biting rates of rural and forest mosquitoes on humans are taken to be β1 = β2 = 0.3. Figure 5 shows that all populations reach the disease- free equilibrium with the disease disappearing from vector populations faster than the human population. For the case R0 > 1, we consider q2 > 0 and q3 > 0 and the biting rate of rural mosquitoes and of forest mosquitoes on humans to be β1 = 0.35, β2 = 0.8, respectively. The unique endemic equilibrium for this case is given by: Z∗ := (0.05867, 0.00054, 0.940789, 0.99808, 0.00192, 0.99890, 0.00109) and the disease dynamics are illustrated in Figure 6. It shows that the solution exhibits oscillations before reaching its steady state. Now, we present numerical simulations for the effects of variations of κ1 and κ2. In Figure 7, we change the Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 11/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika Fig. 5: Phase diagram when R0 < 1 with parameter values taken to be α1 = 2, α2 = 3, κ1 = 0.3, κ2 = 0.2, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 0.07, µU = 0.07, β1 = 0.3, β2 = 0.3, λ = 0.01, γ = 0.16, θ1 = 0.2, θ2 = 0.3. fraction of susceptible humans moving to forest area κ1 and fix the other parameters as listed in Table II. We note that increasing the values of κ1 leads to a slight increase in the maximum of both infected humans and infected vectors in forest areas. The infections reach their maximum and their endemic steady states slightly earlier as the movement of susceptible humans increases. Figure 8 illustrates the effect of varying the propor- tion of infected humans moving to forest area κ2 and fixing all other parameters. It shows that increasing the proportion of infected humans moving to forest areas has the effect of increasing the number of infected vectors in forest areas and the number of infected humans. The time it takes to reach the maximum Fig. 6: Phase diagram when R0 > 1 with parameter values taken to be α1 = 2, α2 = 3, κ1 = 0.3, κ2 = 0.2, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 0.025, µU = 0.025, β1 = 0.35, β2 = 0.8, λ = 0.235, γ = 0.07, θ1 = 0.33, θ2 = 0.3. number of infections remains the same for the vector population and it becomes slightly earlier for the human population as κ2 increases. Note that when κ2 = 0 the number of infected mosquitoes in the forest area reaches zero, which means that the disease will disappear from the forest since the model assumes that infected humans are the only source of infection for the vector population in the forest. However, there have been some reports of ZIKV being found in non-human primates, raising the possibility that they could act as reservoirs [30, 31]. These sources of infection will be considered in future works. V. CONCLUSION A mathematical model of ZIKV disease including hu- man movement and three transmission routes, namely, Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 12/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika 0 50 100 150 t 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 I H (t ) 1 =0.25 1 =0.3 1 =0.35 1 =0.4 0 50 100 150 t 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 I U (t ) 1 =0.25 1 =0.3 1 =0.35 1 =0.4 Fig. 7: Number of infected populations for different values of κ1, where other parameters are fixed κ2 = 0.2, α1 = 2, α2 = 3, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 0.07, µU = 0.07, β1 = 0.35, β2 = 0.9, λ = 0.235, γ = 0.16, θ1 = 0.33, θ2 = 0.3. human-to-human transmission, vector transmission, and vertical transmission has been proposed. The model has been analyzed and studied to investigate the role of human movement from rural areas to forest areas on the spread of ZIKV. The positivity of the solution and the boundedness of the invariant region were discussed. The basic reproduction number R0 was computed and expressed in terms of reproduction numbers related to the interactions between humans RHH, between human and vector in rural area RHV and between human and vector in forest area RHU . It was found that the threshold of the disease which occurs at R0 = 1 is equivalent to RHH + RHV + RHU = 1 and hence all transmission routes need to be controlled to reduce the spread of the disease. 0 50 100 150 200 250 300 t 0 0.05 0.1 0.15 I H (t ) 2 =0 2 =0.05 2 =0.1 2 =0.15 0 50 100 150 200 250 300 t 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 I U (t ) 2 =0 2 =0.05 2 =0.1 2 =0.15 Fig. 8: Number of infected populations for different values of κ2, where other parameters are fixed κ1 = 0.3, α1 = 2, α2 = 3, ε1 = 0.67, ε2 = 0.06, ε3 = 0.06, µH = 0.00004, µV = 0.07, µU = 0.07, β1 = 0.35, β2 = 0.9, λ = 0.235, γ = 0.16, θ1 = 0.33, θ2 = 0.3. Sensitivity analysis of R0 was carried out and it showed that R0 is sensitive to almost all model pa- rameters either positively or negatively, except the pa- rameters κ1 and κ2, representing the proportions of susceptible and infected humans moving to forest areas, respectively, where their signs of sensitivity indices were found to depend on the biting rates when fixing all other parameters. This dependence has been calculated and illustrated graphically. It has been found that the indices are both positive when the forest mosquito’s biting rate is high and the rural mosquito’s biting rate is small up to a certain limit. However, this positive effect on R0 was found to be very small, with κ2 having a higher effect than κ1. The most positive influential parameters are the biting rate Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 13/14 https://doi.org/10.55630/j.biomath.2022.12.149 Al-Maqrashi et al, Investigating the Role of Mobility between Rural Areas and Forests on the Spread of Zika of rural and forest mosquitoes on humans, while the recovery rate of humans has the most negative impact. Then, the local and global stability of the disease- free equilibrium was derived whenever R0 is less than unity. Furthermore, the system was shown to possess a unique endemic equilibrium under certain conditions, and it is locally asymptotically stable when R0 is greater than unity since the direction of the bifurcation was found to be forward. The bifurcation analysis was presented both analytically and graphically. Finally, numerical simulations were presented to demonstrate the obtained theoretical results. They confirmed that the human movement from rural areas to forests has a small effect on increasing the infected human and vector populations, with the movement of infected humans having a higher effect than the movement of susceptible humans. REFERENCES [1] G. S. Campos, A. C. Bandeira, S. I. Sardi, “Zika Virus Outbreak, Bahia, Brazil”, Emerging Infectious Diseases, 21(10):1885- 1886, 2015. [2] V.-M. Cao-Lormeau, A. Blake, S. Mons, S. Lastère, C. Roche, J. Vanhomwegen, et al., “Guillain-Barré Syndrome outbreak associated with Zika virus infection in French Polynesia: a case- control study”, The Lancet, 387:1531-1539, 2016. [3] N. K. Goswami, A. K. Srivastav, M. Ghosh, B. Shanmukha, “Mathematical modeling of zika virus disease with nonlinear incidence and optimal control”, Journal of Physics: Conference Series, 1000:012114, 2018. [4] S. K. Biswas, U. Ghosh, S. Sarkar, “Mathematical model of zika virus dynamics with vector control and sensitivity analysis”, Infectious Disease Modelling, 5:23-41, 2020. [5] A. Troncoso, “Zika threatens to become a huge worldwide pandemic”, Asian Pacific Journal of Tropical Biomedicine, 6(6):520-527, 2016. [6] L. Schuler-Faccini, E. M. Ribeiro, I. M. L. Feitosa, D. D. G. Horovitz, D. P. Cavalcanti, A. Pessoa, et al., “Possible Association Between Zika Virus Infection and Microcephaly – Brazil, 2015”, MMWR Morbidity and Mortality Weekly Report, 65:59-62, 2016. [7] V. Heang, C. Y. Yasuda, L. Sovann, A. D. Haddow, A. Travassos da Rosa, R. B. Tesh, M. R. Kasper, “Zika Virus Infection, Cambodia, 2010”, Emerging Infectious Diseases, 18(2):349- 351, 2012. [8] E. Bonyah, M. A. Khan, K. O. Okosun, S. Islam, “A theoretical model for Zika virus transmission”, PLoS ONE, 12(10):e0185540, 2017. [9] F. B. Agusto, S. Bewick, W. F. Fagan, “Mathematical model of Zika virus with vertical transmission”, Infectious Disease Modelling, 2(2):244-267, 2017. [10] Z. C. Hlaing, C. Kamiyama, O. Saito, “Interaction between Rural People’s Basic Needs and Forest Products: A Case Study of the Katha District of Myanmar”, International Journal of Forestry Research, 2017:2105012, 2017. [11] Centers for Disease Control and Prevention (CDC), Zika Virus – Prevention and Transmission, https://www.cdc.gov/zika/prevention/, [16/12/2022]. [12] P. Brasil, Z. Vasconcelos, T. Kerin, C. R. Gabaglia, I. P. Ribeiro, M. C. Bonaldo, et al., “Zika virus vertical transmission in children with confirmed antenatal exposure”, Nature Commu- nications, 11:3510, 2020. [13] S. Thangamani, J. Huang, C. E. Hart, H. Guzman, R. B. Tesh, “Vertical Transmission of Zika Virus in Aedes aegypti Mosquitoes”, The American Journal of Tropical Medicine and Hygiene, 95(5):1169-1173, 2016. [14] E. O. Alzahrani, W. Ahmad, M. A. Khan, S. J. Malebary, “Optimal Control Strategies of Zika Virus Model with Mutant”, Communications in Nonlinear Science and Numerical Simula- tion, 93:105532, 2021. [15] Pan American Health Organization (PAHO), Zika Virus, https://www.paho.org/en/topics/zika, [16/12/2022]. [16] D. P. Shutt, C. A. Manore, S. Pankavich, A. T. Porter, S. Y. Del Valle, “Estimating the reproductive number, total outbreak size, and reporting rates for Zika epidemics in South and Central America”, Epidemics, 21:63-79, 2017. [17] O. Maxian, A. Neufeld, E. J. Talis, L. M. Childs, J. C. Black- wood, “Zika virus dynamics: When does sexual transmission matter?”, Epidemics, 21:48-55, 2017. [18] A. Dénes, M. A. Ibrahim, L. Oluoch, M. Tekeli, T. Tekeli, “Impact of weather seasonality and sexual transmission on the spread of Zika fever”, Scientific Reports, 9:17055, 2019. [19] Z. Lai, T. Zhou, J. Zhou, S. Liu, Y. Xu, J. Gu, G. Yan, X.-G. Chen, “Vertical transmission of zika virus in Aedes albopictus”, PLoS Neglected Tropical Diseases, 14(10):e0008776, 2020. [20] P. van den Driessche, J. Watmough, “Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission”, Mathematical Biosciences, 180(1-2):29- 48, 2002. [21] H. S. Rodrigues, M. T. T. Monteiro, D. F. M. Torres, “Sensitivity Analysis in a Dengue Epidemiological Model”, Conference Papers in Science, 2013:721406, 2013. [22] J. Martinez, A. Showering, C. Oke, R. T. Jones, J. G. Logan, “Differential attraction in mosquito-human interactions and im- plications for disease control”, Philosophical Transactions of the Royal Society B, 376(1818), 2021. [23] D. A. M. Villela, L. S. Bastos, L. M. de Carvalho, O. G. Cruz, M. F. C. Gomes, B. Durovni, M. C. Lemos, V. Saraceni, F. C. Coelho, C. T. Codeço, “Zika in Rio de Janeiro: Assessment of basic reproduction number and comparison with dengue out- breaks”, Epidemiology and Infection, 145(8):1649-1657, 2017. [24] D. K. Sambariya, R. Prasad, “Routh Stability Array Method based reduced model of Single Machine Infinite Bus with Power System Stabilizer”, International Conference on Emerg- ing Trends in Electrical, Communication and Information Tech- nologies (ICECIT-2012), 2012. [25] C. Castillo-Chavez, Z. Feng, W. Huang, “On the computation of R0 and its role on global stability”, Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Intro- duction, 1:229-250, 2002, download link. [26] D. R. Curtiss, “Recent Extentions of Descartes’ Rule of Signs”, Annals of Mathematics, Second Series, 19(4):251-278, 1918. [27] J. H. Hubbard, B. H. West, Differential Equations: A Dynam- ical Systems Approach. Higher Dimensional Systems, Texts in Applied Mathematics, 18, Springer-Verlag, New York, 1995. [28] A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G.E. Meijer, B. Sautois, “New features of the software MatCont for bifurcation analysis of dynamical systems”, Mathematical and Computer Modelling of Dynamical Systems, 14(2):147-175, 2008. [29] C. Castillo-Chavez, B. Song, “Dynamical Models of Tubercu- losis and Their Applications”, Mathematical Biosciences and Engineering, 1(2):361-404, 2004. [30] D. M. Dudley, M. T. Aliota, E. Mohr, A. M. Weiler, G. Lehrer- Brey, K. L. Weisgrau, et al., “A rhesus macaque model of Asia lineage Zika virus infection”, bioRxiv 046334, 2016. [31] N. N. Haese, V. H. J. Roberts, A. Chen, D. N. Streblow, T. K. Morgan, A. J. Hirsch, “Nonhuman Primate Models of Zika Virus Infection and Disease during Pregnancy”, Viruses, 13(10):2088, 2021. Biomath 11 (2022), 2212149, https://doi.org/10.55630/j.biomath.2022.12.149 14/14 https://doi.org/10.3201/eid2110.150847 https://doi.org/10.3201/eid2110.150847 https://doi.org/10.3201/eid2110.150847 https://doi.org/10.1016/S0140-6736(16)00562-6 https://doi.org/10.1016/S0140-6736(16)00562-6 https://doi.org/10.1016/S0140-6736(16)00562-6 https://doi.org/10.1016/S0140-6736(16)00562-6 https://doi.org/10.1088/1742-6596/1000/1/012114 https://doi.org/10.1088/1742-6596/1000/1/012114 https://doi.org/10.1088/1742-6596/1000/1/012114 https://doi.org/10.1088/1742-6596/1000/1/012114 https://doi.org/10.1016/j.idm.2019.12.001 https://doi.org/10.1016/j.idm.2019.12.001 https://doi.org/10.1016/j.idm.2019.12.001 https://doi.org/10.1016/j.apjtb.2016.04.004 https://doi.org/10.1016/j.apjtb.2016.04.004 https://doi.org/10.1016/j.apjtb.2016.04.004 https://doi.org/10.15585/mmwr.mm6503e2 https://doi.org/10.15585/mmwr.mm6503e2 https://doi.org/10.15585/mmwr.mm6503e2 https://doi.org/10.15585/mmwr.mm6503e2 https://doi.org/10.15585/mmwr.mm6503e2 https://doi.org/10.3201/eid1802.111224 https://doi.org/10.3201/eid1802.111224 https://doi.org/10.3201/eid1802.111224 https://doi.org/10.3201/eid1802.111224 https://doi.org/10.1371/journal.pone.0185540 https://doi.org/10.1371/journal.pone.0185540 https://doi.org/10.1371/journal.pone.0185540 https://doi.org/10.1016/j.idm.2017.05.003 https://doi.org/10.1016/j.idm.2017.05.003 https://doi.org/10.1016/j.idm.2017.05.003 https://doi.org/10.1155/2017/2105012 https://doi.org/10.1155/2017/2105012 https://doi.org/10.1155/2017/2105012 https://doi.org/10.1155/2017/2105012 https://www.cdc.gov/zika/prevention/ https://www.cdc.gov/zika/prevention/ https://www.cdc.gov/zika/prevention/ https://doi.org/10.1038/s41467-020-17331-0 https://doi.org/10.1038/s41467-020-17331-0 https://doi.org/10.1038/s41467-020-17331-0 https://doi.org/10.1038/s41467-020-17331-0 https://doi.org/10.4269/ajtmh.16-0448 https://doi.org/10.4269/ajtmh.16-0448 https://doi.org/10.4269/ajtmh.16-0448 https://doi.org/10.4269/ajtmh.16-0448 https://doi.org/10.1016/j.cnsns.2020.105532 https://doi.org/10.1016/j.cnsns.2020.105532 https://doi.org/10.1016/j.cnsns.2020.105532 https://doi.org/10.1016/j.cnsns.2020.105532 https://www.paho.org/en/topics/zika https://www.paho.org/en/topics/zika https://doi.org/10.1016/j.epidem.2017.06.005 https://doi.org/10.1016/j.epidem.2017.06.005 https://doi.org/10.1016/j.epidem.2017.06.005 https://doi.org/10.1016/j.epidem.2017.06.005 https://doi.org/10.1016/j.epidem.2017.06.003 https://doi.org/10.1016/j.epidem.2017.06.003 https://doi.org/10.1016/j.epidem.2017.06.003 https://doi.org/10.1038/s41598-019-53062-z https://doi.org/10.1038/s41598-019-53062-z https://doi.org/10.1038/s41598-019-53062-z https://doi.org/10.1371/journal.pntd.0008776 https://doi.org/10.1371/journal.pntd.0008776 https://doi.org/10.1371/journal.pntd.0008776 https://doi.org/10.1016/S0025-5564(02)00108-6 https://doi.org/10.1016/S0025-5564(02)00108-6 https://doi.org/10.1016/S0025-5564(02)00108-6 https://doi.org/10.1016/S0025-5564(02)00108-6 https://doi.org/10.1155/2013/721406 https://doi.org/10.1155/2013/721406 https://doi.org/10.1155/2013/721406 https://doi.org/10.1098/rstb.2019.0811 https://doi.org/10.1098/rstb.2019.0811 https://doi.org/10.1098/rstb.2019.0811 https://doi.org/10.1098/rstb.2019.0811 https://doi.org/10.1017/S0950268817000358 https://doi.org/10.1017/S0950268817000358 https://doi.org/10.1017/S0950268817000358 https://doi.org/10.1017/S0950268817000358 https://doi.org/10.1017/S0950268817000358 https://doi.org/10.13140/RG.2.1.4041.8325 https://doi.org/10.13140/RG.2.1.4041.8325 https://doi.org/10.13140/RG.2.1.4041.8325 https://doi.org/10.13140/RG.2.1.4041.8325 https://doi.org/10.13140/RG.2.1.4041.8325 https://books.google.com/books?id=pR4CqiTSTMwC https://books.google.com/books?id=pR4CqiTSTMwC https://www.researchgate.net/profile/Carlos-Castillo-Chavez/publication/228915276_ON_THE_COMPUTATION_OF_R_AND_ITS_ROLE_ON_GLOBAL_STABILITY/links/0fcfd50a569d32c708000000/ON-THE-COMPUTATION-OF-R-AND-ITS-ROLE-ON-GLOBAL-STABILITY.pdf https://doi.org/10.2307/1967494 https://doi.org/10.2307/1967494 https://doi.org/10.1007/978-1-4612-4192-8 https://doi.org/10.1007/978-1-4612-4192-8 https://doi.org/10.1007/978-1-4612-4192-8 https://doi.org/10.1080/13873950701742754 https://doi.org/10.1080/13873950701742754 https://doi.org/10.1080/13873950701742754 https://doi.org/10.1080/13873950701742754 https://doi.org/10.3934/mbe.2004.1.361 https://doi.org/10.3934/mbe.2004.1.361 https://doi.org/10.3934/mbe.2004.1.361 https://doi.org/10.1101/046334 https://doi.org/10.1101/046334 https://doi.org/10.1101/046334 https://doi.org/10.3390/v13102088 https://doi.org/10.3390/v13102088 https://doi.org/10.3390/v13102088 https://doi.org/10.3390/v13102088 https://doi.org/10.55630/j.biomath.2022.12.149 Introduction Model Description Model Analysis Positivity of Solutions and Positively Invariant Set The Basic Reproduction Number Sensitivity Analysis of the Basic Reproduction Number Local Stability of the DFE Global Stability of the DFE Existence of Endemic Equilibrium Bifurcation Analysis Numerical Analysis Conclusion References