# Basic Properties and Qualitative Dynamics of a Vector-Borne Disease Model with Vector Stages and Vertical Transmission

Authors: Sansao A. Pedro
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

This work systematically discusses basic properties and qualitative dynamics of vector-borne disease models, particularly those with vertical transmission in the vector population. Examples of disease include Dengue and Rift Valley fever which are endemic in Sub-Saharan Africa, and understanding of the dynamics underlying their transmission is central for providing critical informative indicators useful for guiding control strategies. Of particular interest is the applicability and derivation of relevant population and epidemic thresholds and their relationships with vertical infection. This study demonstrates how the failure ofR0 derived using the next-generation method compounds itself when varying vertical transmission efficiency, and it shows that the host type reproductive number gives the correct R0. Further, novel relationships between the host type reproductive number, vertical infection, and ratio of female mosquitoes to host are established and discussed. Analytical results of the model with vector stages show that the quantities Q0, Q0v, and R0c, which represent the vector colonization threshold, the average number of female mosquitoes produced by a single infected mosquito, and effective reproductive number, respectively, provide threshold conditions that determine the establishment of the vector population and invasion of the disease. Numerical simulations are also conducted to confirm and extend the analytical results. The findings imply that while vertical infection increases the size of an epidemic, it reduces its duration, and control efforts aimed at reducing the critical thresholds Q0, Q0v, and R0c to below unity are viable control strategies.


---

## Body

## 1. Introduction

Vector-borne diseases have been the scourge of man and animals since the beginning of time [1]. Today, vector-borne diseases account for over 17% of all infectious diseases causing more than 1 million deaths annually, and their distribution is mainly determined by a complex dynamic of environmental and social factors [2]. In spite of all these inherent complexities, mathematical models have been used to translate assumptions concerning biological, environmental, and social aspects into mathematical structures, linking biological processes of transmission and dynamics of infection at population level. Such dynamic models have impacted both our understanding of epidemic spread and public health planning (for more details see [3–5] and references therein).In this study our particular interest is in investigating qualitative properties of epidemic models of mosquito-borne diseases in which the vector is of generaAedes. These mosquito species are known to transmit many vector-borne diseases of vast epidemiological importance including Dengue fever and Rift Valley fever (RVF), just to mention a few. These diseases are endemic in Sub-Saharan Africa with pronounced health and economic impacts on domestic animals and humans. In fact, losses due to RVF can reach millions of dollars during periods of disease outbreaks [6, 7]. An interesting phenomenon underlying many vector-borne diseases is their ability to persist year round, fluctuating seasonally but only falling to zero during some months. Hence, an important question is, how does the virus maintain itself in nature? For RVF it has been hypothesized that RVF virus (RVFV) is maintained through transovarial transmission inAedes mosquito eggs [8].Aedes eggs need be dry for several days before they can mature. After maturing, they hatch during the next flooding event large enough to cover them with water [9–11]. The eggs have high desiccation resistance and can survive dry conditions in a dormant form for months to years [12–15]. Thus, the emergence of adultAedes mosquitoes from infected eggs can reintroduce RVF in livestock at the beginning of the rainy season, before other mosquitoes species amplify it further [16]. For instance, in eastern and southern Africa there is more and more evidence of disease activities between outbreaks [7, 17–20] highlighting the role of vertical transmission for initial disease spread and endemicity.In epidemiology disease spread and persistence are measured through quantities known as epidemic thresholds. Their derivation and characterization are one of the most important results of mathematical epidemic models. The basic reproductive number,R0, is the most critical epidemic threshold given its applicability and suitability for deciding whether an outbreak will occur or fade out, making it essential for guiding disease control efforts. However, the derivation of this epidemic threshold in vector-borne disease models in particular suffers from a lack of uniqueness and it fails to give the correct average number of expected secondary infections produced by one infected individual [21]. This failure is more likely to compound itself when vertical transmission mode is included in the transmission model, since the resulting R0 comes as the sum of the vertical and horizontal transmission components, if the next-generation method is used. Previous mathematical models have made a significant attempt in including vertical infection in modelling vector-borne diseases [22–28], but none of them discuss how the failure of R0 compounds itself in the presence of vertical transmission. Therefore, the present work aims to discuss some relevant basic properties of vector-borne disease models when vertical infection is taken into account and their implications for disease control efforts. Further, our goal is to derive new epidemic thresholds useful for guiding control efforts in the settings of vector-borne disease models with vector stages that include vertical transmission mode.We formulate two models, one simple but realistic and the other more complex with vector stages. The first is an extension of the one proposed by Ross [29] and popularized by Macdonald [30] and Anderson [31]. The model is used to discuss system properties such as the asymmetric relationship between the host-to-vector and vector-to-host reproductive numbers. In addition, we highlight how to derive epidemic thresholds useful for guiding disease control efforts and discuss their relationships with vertical transmission efficiency. It is shown that the model has two model equilibria, namely, the disease-free and the endemic, and Lyapunov function theory is used to establish their global qualitative dynamics. The second is an extension of the basic model, where the dynamics of both aquatic and adult mosquitoes are modelled explicitly. In this model we let the populations of aquatic and adults vary with time but be limited by their respective carrying capacity. The inclusion of the explicit vector submodel allows for derivation of critical thresholds such as the reproductive number for both the vector population and the disease system. Then, these thresholds are used to determine global qualitative dynamics of both the disease-free and endemic equilibria.The paper is set out as follows. In Section2 we formulate and discuss the two model systems. In Section 3 we provide the epidemic threshold theorems regarding both the vector population and disease equilibria. We also discuss important model properties and how to derive and identify model epidemic thresholds useful for guiding disease control efforts. Furthermore, numerical simulations are carried out to investigate the influence of the key parameters on the spread of the disease (taking RVF as disease example), to support analytical analyses and conclusions and illustrate possible behavioural scenarios of the model with vector stages. Finally, in Section 4 we present a short discussion of the results and their biological implications.


## 2. Materials and Methods

For human and animal diseases, horizontal transmission typically occurs through direct or indirect physical contact with infectious hosts, or through disease vectors such as mosquitoes, ticks, or other biting insects. Among mosquito vectors vertical transmission is often through eggs. Of particular interest are female mosquitoes of generaAedes which transmit the virus to their eggs. These eggs have some adaptive behaviour which allows them to stay dormant in nature for relatively long periods. Although vertical transmission also occurs among vectors involved in the transmission of Dengue disease, RVF is the disease for which the model is a good approximation. In particular, parameter values related to RVF are used to illustrate the dynamics of the disease numerically.


### 2.1. Host-Vector Basic Model with Vertical Transmission

LetNh=Sh+Ih+Rh and Nv=Sv+Iv denote the total host and vector populations sizes, respectively. We assume that individuals at each compartment mix homogeneously and each mosquito bites each individual host at a constant rate α/Nh, where α is the biting rate per unit time. Let phv be a probability of successful infection transmission from an infected mosquito to a susceptible host and pvh be a probability of successful infection transmission from an infected host to a susceptible mosquito per bite. Thus, the forces of infection are as follows, λhv=phv(α/Nh)Iv and λvh=phv(α/Nh)Ih. Hosts are recruited into the population at per capita rate μh which is proportional to the total population and leave each compartment through death. Noninfected mosquitoes join the susceptible compartment at rate μv(Nv-qvIv) while vertically infected mosquitoes join the infected class at rate μvqvIv. Assuming constant population sizes, that is, births equal to deaths, Rh can be obtained when both Sh and Ih are known. In the same way, Sv can be obtained when Iv is known. Thus, the expressions for both Rh and Sv can be omitted and the system can be written in terms of proportions:(1)dShdt=μh1-Sh-phvαmIvSh,dIhdt=phvαmIvSh-γh+μhIh,dIvdt=pvhαIh1-Iv-μv1-qvIv,where qv∈(0,1) is the probability of vertical infection, γh the rate at which infected hosts recover from infection, and m denotes the ratio of female mosquitoes to hosts.


### 2.2. Model with Vector Stages and Vertical Transmission

Here we extend the basic model to include vector stages. Partial results of the resulting model without vertical infection have been obtained in [32]. Our aim is to extend their model and analysis by investigating global dynamics of all model equilibria and examine the extent to which vertical infection alters the dynamics of the system. The mosquito population is divided into aquatic (eggs, larvae, and pupae) and terrestrial (adults) subpopulations with K1 and K2 being their carrying capacity, respectively. The parameter K1 represents the larval maximal capacity limited by the availability of breeding sites while K2 is the maximal capacity of adult mosquitoes limited by factors conditioning their survival such as high altitudes and high temperatures. Further, the aquatic subpopulation is divided into epidemiological classes, susceptible (Ls), and infected (Li) while adults are divided also into susceptible (Sv) and infected (Iv). The per capita oviposition rate is r(1-L/K1), where r is the intrinsic oviposition rate and L=Ls+Li. Aquatic mosquitoes emerge as adults at a per capita rate θ where the proportion 1-qv emerge noninfected while the remainder are infected. Disease transmission dynamics between vector and host populations remain the same as in the basic model. As a result the following nondimensional system of ordinary differential equations represents the model that governs the temporal evolution of the disease:(2)dLsdt=δr1-LNv-qvIv-μL+θLs,dLidt=δr1-LqvIv-μL+θLi,dSvdt=θδLs1-Nv-λvhIhSv-μvSv,dIvdt=θδLi1-Nv+λvhIhSv-μvIv,dShdt=μh1-Sh-λhvIvSh,dIhdt=λhvIvSh-γh+μhIh, where L=Ls+Li,δ=K2/K1, λvh=pvhα, λhv=phvαm, m=K2/Nh, Sv=Sv/K2, and Iv=Iv/K2.LetLs(0),Li(0),Sv(0),Iv(0),Sh(0),Ih(0)⩾0∈R+6 be the initial conditions of system (2). It is easy to check that the feasible region for (2) is the positive orthant of R6 and that the closed set(3)Φ=Ls,Li,Sv,Iv,Sh,Ih∈R+6:Ls+Li⩽1,Sv+Iv⩽1,Sh+Ih⩽1is positively invariant for system (2).


#### 2.2.1. Positivity of Solutions

Lemma 1 (see [33]).
Let us denoteu=(ui)i=1,n¯ and consider the function f:R+×Rn→Rn continuous with respect to t,u and Lipschitz with respect to u. If f(t,u)≥0 for (t,u)∈R+×Rn, with ui=0, then, for every u0∈R+n, there exists T>0 such that the solution to (4)dudt=ft,u,ut0=u0exists and is unique and positive with value in R+n and defined on some interval [0,T). If T<∞, then (5)limsupt→T∑i=1nui=+∞.Theorem 2.
The solution set{Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)} of the model (2) exists and is unique and positive for t>0.Proof.
LetX=(Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)) and denote the function f=(fi)i=1,6¯ such that(6)f1Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LNv-qvIv-μL+θLs,f2Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LqvIv-μL+θLi,f3Lst,Lit,Svt,Ivt,Sht,Iht=θδLs1-Nv-λvhIhSv-μvSv,f4Lst,Lit,Svt,Ivt,Sht,Iht=θδLi1-Nv+λvhIhSv-μvIv,f5Lst,Lit,Svt,Ivt,Sht,Iht=μh1-Sh-λhvIvSh,f6Lst,Lit,Svt,Ivt,Sht,Iht=λhvIvSh-γh+μhIh.
Since the functionf is continuous and Lipschitz continuous with respect to X, according to Picard’s theorem, there exists T0>0 such that the solution to (2) exists locally at least on an interval of this form [0,T0]. Further, considering the initial condition X1=X(T0) at t0=T0 and using Picard’s theorem, it follows that there exists T0≤T1∈R+ such that the solution to (2) exists and is unique on [T0,T1]. Since f is continuous and differentiable, the solution of (2) with a given initial condition is unique. Therefore, the solutions of (2) obtained on [0,T0] and on [T0,T1] form the unique solution of (2) on IT1=[0,T1] with the initial condition X0 at t0=0. Repeating this process again and again, we end up with the maximal forward interval of existence for the solutions of (2), say IT^=[0,T^) with T^>0. Furthermore, for X∈R+6, (7)fiLst,Lit,Svt,Ivt,Sht,0≥0. Therefore, the solutions of (2) on IT^ are positive, according to Lemma 1.
Finally, according to Theorem2, the solutions to (2) are bounded on [0,T^). In other words, they do not blow up on any finite interval of R+. It follows that, according to Lemma 1, the solution of (2) exists for all time. Hence for any initial condition in R+6, system (2) possesses a unique and positive solution in R+6.


## 2.1. Host-Vector Basic Model with Vertical Transmission

LetNh=Sh+Ih+Rh and Nv=Sv+Iv denote the total host and vector populations sizes, respectively. We assume that individuals at each compartment mix homogeneously and each mosquito bites each individual host at a constant rate α/Nh, where α is the biting rate per unit time. Let phv be a probability of successful infection transmission from an infected mosquito to a susceptible host and pvh be a probability of successful infection transmission from an infected host to a susceptible mosquito per bite. Thus, the forces of infection are as follows, λhv=phv(α/Nh)Iv and λvh=phv(α/Nh)Ih. Hosts are recruited into the population at per capita rate μh which is proportional to the total population and leave each compartment through death. Noninfected mosquitoes join the susceptible compartment at rate μv(Nv-qvIv) while vertically infected mosquitoes join the infected class at rate μvqvIv. Assuming constant population sizes, that is, births equal to deaths, Rh can be obtained when both Sh and Ih are known. In the same way, Sv can be obtained when Iv is known. Thus, the expressions for both Rh and Sv can be omitted and the system can be written in terms of proportions:(1)dShdt=μh1-Sh-phvαmIvSh,dIhdt=phvαmIvSh-γh+μhIh,dIvdt=pvhαIh1-Iv-μv1-qvIv,where qv∈(0,1) is the probability of vertical infection, γh the rate at which infected hosts recover from infection, and m denotes the ratio of female mosquitoes to hosts.


## 2.2. Model with Vector Stages and Vertical Transmission

Here we extend the basic model to include vector stages. Partial results of the resulting model without vertical infection have been obtained in [32]. Our aim is to extend their model and analysis by investigating global dynamics of all model equilibria and examine the extent to which vertical infection alters the dynamics of the system. The mosquito population is divided into aquatic (eggs, larvae, and pupae) and terrestrial (adults) subpopulations with K1 and K2 being their carrying capacity, respectively. The parameter K1 represents the larval maximal capacity limited by the availability of breeding sites while K2 is the maximal capacity of adult mosquitoes limited by factors conditioning their survival such as high altitudes and high temperatures. Further, the aquatic subpopulation is divided into epidemiological classes, susceptible (Ls), and infected (Li) while adults are divided also into susceptible (Sv) and infected (Iv). The per capita oviposition rate is r(1-L/K1), where r is the intrinsic oviposition rate and L=Ls+Li. Aquatic mosquitoes emerge as adults at a per capita rate θ where the proportion 1-qv emerge noninfected while the remainder are infected. Disease transmission dynamics between vector and host populations remain the same as in the basic model. As a result the following nondimensional system of ordinary differential equations represents the model that governs the temporal evolution of the disease:(2)dLsdt=δr1-LNv-qvIv-μL+θLs,dLidt=δr1-LqvIv-μL+θLi,dSvdt=θδLs1-Nv-λvhIhSv-μvSv,dIvdt=θδLi1-Nv+λvhIhSv-μvIv,dShdt=μh1-Sh-λhvIvSh,dIhdt=λhvIvSh-γh+μhIh, where L=Ls+Li,δ=K2/K1, λvh=pvhα, λhv=phvαm, m=K2/Nh, Sv=Sv/K2, and Iv=Iv/K2.LetLs(0),Li(0),Sv(0),Iv(0),Sh(0),Ih(0)⩾0∈R+6 be the initial conditions of system (2). It is easy to check that the feasible region for (2) is the positive orthant of R6 and that the closed set(3)Φ=Ls,Li,Sv,Iv,Sh,Ih∈R+6:Ls+Li⩽1,Sv+Iv⩽1,Sh+Ih⩽1is positively invariant for system (2).


### 2.2.1. Positivity of Solutions

Lemma 1 (see [33]).
Let us denoteu=(ui)i=1,n¯ and consider the function f:R+×Rn→Rn continuous with respect to t,u and Lipschitz with respect to u. If f(t,u)≥0 for (t,u)∈R+×Rn, with ui=0, then, for every u0∈R+n, there exists T>0 such that the solution to (4)dudt=ft,u,ut0=u0exists and is unique and positive with value in R+n and defined on some interval [0,T). If T<∞, then (5)limsupt→T∑i=1nui=+∞.Theorem 2.
The solution set{Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)} of the model (2) exists and is unique and positive for t>0.Proof.
LetX=(Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)) and denote the function f=(fi)i=1,6¯ such that(6)f1Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LNv-qvIv-μL+θLs,f2Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LqvIv-μL+θLi,f3Lst,Lit,Svt,Ivt,Sht,Iht=θδLs1-Nv-λvhIhSv-μvSv,f4Lst,Lit,Svt,Ivt,Sht,Iht=θδLi1-Nv+λvhIhSv-μvIv,f5Lst,Lit,Svt,Ivt,Sht,Iht=μh1-Sh-λhvIvSh,f6Lst,Lit,Svt,Ivt,Sht,Iht=λhvIvSh-γh+μhIh.
Since the functionf is continuous and Lipschitz continuous with respect to X, according to Picard’s theorem, there exists T0>0 such that the solution to (2) exists locally at least on an interval of this form [0,T0]. Further, considering the initial condition X1=X(T0) at t0=T0 and using Picard’s theorem, it follows that there exists T0≤T1∈R+ such that the solution to (2) exists and is unique on [T0,T1]. Since f is continuous and differentiable, the solution of (2) with a given initial condition is unique. Therefore, the solutions of (2) obtained on [0,T0] and on [T0,T1] form the unique solution of (2) on IT1=[0,T1] with the initial condition X0 at t0=0. Repeating this process again and again, we end up with the maximal forward interval of existence for the solutions of (2), say IT^=[0,T^) with T^>0. Furthermore, for X∈R+6, (7)fiLst,Lit,Svt,Ivt,Sht,0≥0. Therefore, the solutions of (2) on IT^ are positive, according to Lemma 1.
Finally, according to Theorem2, the solutions to (2) are bounded on [0,T^). In other words, they do not blow up on any finite interval of R+. It follows that, according to Lemma 1, the solution of (2) exists for all time. Hence for any initial condition in R+6, system (2) possesses a unique and positive solution in R+6.


## 2.2.1. Positivity of Solutions

Lemma 1 (see [33]).
Let us denoteu=(ui)i=1,n¯ and consider the function f:R+×Rn→Rn continuous with respect to t,u and Lipschitz with respect to u. If f(t,u)≥0 for (t,u)∈R+×Rn, with ui=0, then, for every u0∈R+n, there exists T>0 such that the solution to (4)dudt=ft,u,ut0=u0exists and is unique and positive with value in R+n and defined on some interval [0,T). If T<∞, then (5)limsupt→T∑i=1nui=+∞.Theorem 2.
The solution set{Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)} of the model (2) exists and is unique and positive for t>0.Proof.
LetX=(Ls(t),Li(t),Sv(t),Iv(t),Sh(t),Ih(t)) and denote the function f=(fi)i=1,6¯ such that(6)f1Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LNv-qvIv-μL+θLs,f2Lst,Lit,Svt,Ivt,Sht,Iht=δr1-LqvIv-μL+θLi,f3Lst,Lit,Svt,Ivt,Sht,Iht=θδLs1-Nv-λvhIhSv-μvSv,f4Lst,Lit,Svt,Ivt,Sht,Iht=θδLi1-Nv+λvhIhSv-μvIv,f5Lst,Lit,Svt,Ivt,Sht,Iht=μh1-Sh-λhvIvSh,f6Lst,Lit,Svt,Ivt,Sht,Iht=λhvIvSh-γh+μhIh.
Since the functionf is continuous and Lipschitz continuous with respect to X, according to Picard’s theorem, there exists T0>0 such that the solution to (2) exists locally at least on an interval of this form [0,T0]. Further, considering the initial condition X1=X(T0) at t0=T0 and using Picard’s theorem, it follows that there exists T0≤T1∈R+ such that the solution to (2) exists and is unique on [T0,T1]. Since f is continuous and differentiable, the solution of (2) with a given initial condition is unique. Therefore, the solutions of (2) obtained on [0,T0] and on [T0,T1] form the unique solution of (2) on IT1=[0,T1] with the initial condition X0 at t0=0. Repeating this process again and again, we end up with the maximal forward interval of existence for the solutions of (2), say IT^=[0,T^) with T^>0. Furthermore, for X∈R+6, (7)fiLst,Lit,Svt,Ivt,Sht,0≥0. Therefore, the solutions of (2) on IT^ are positive, according to Lemma 1.
Finally, according to Theorem2, the solutions to (2) are bounded on [0,T^). In other words, they do not blow up on any finite interval of R+. It follows that, according to Lemma 1, the solution of (2) exists for all time. Hence for any initial condition in R+6, system (2) possesses a unique and positive solution in R+6.


## 3. Results

### 3.1. Analysis of the Basic Host-Vector Model

#### 3.1.1. Model Equilibria and Stability Analysis

The basic host-vector model with vertical transmission exhibits two equilibria, namely, the disease-freeE10 and the endemic E1∗, respectively. At the disease-free equilibrium,(8)E10=Sh0,Ih0,Iv0=1,0,0,both vector and host populations persist but with no disease. The prevalence of the disease is denoted by(9)E1∗=Sh∗,Ih∗,Iv∗=a+μhR0ea+μhR0h,μhR0e-1cR0e+bμh/V,μhR0e-1a+μhR0e,where a=phvαm, b=pvhα, c=γh+μh, V=μv(1-qv), and(10)R0e=abVc=phvpvhα2mμv1-qvγh+μh.The Jacobian matrix of system (1) at E10 is given by(11)J1′=-μh0-a0-ca0b-V,such that the characteristic polynomial of matrix (11) is then given as(12)P1λ=λ3+a2λ2+a1λ+a0,with a2=V+c+μh, a1=Vc(1-R0e)+(V+c)μh, and a0=μhVc(1-R0e). The coefficient a2>0 and both a1,a0 are nonnegative if and only if R0e<1. Hence, all Routh stability criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Furthermore, for R0e=1, (9) becomes the disease-free equilibrium. Therefore, the following result holds.Theorem 3.
The disease-free equilibriumE10=(1,0,0) exists and it is globally asymptotically stable if R0e⩽1.Alternatively, the global stability of the disease-free equilibriumE10 can be established using the following Lyapunov function:(13)VSh,Ih,Iv=u1Sh-Sh0-Sh0lnShSh0+u2Ih+u3Ivwhere u1,u2,u3 are some positive constants. Calculating the derivative of V along the solutions of system (1), we obtain(14)V′=u1Sh-Sh0ShμhSh0-Sh-aShIv+u2aShIv-γh+μhIh+u3bIh1-Iv-μv1-qvIv,(15)=-u1μhSh0-Sh2Sh-u3μv1-qv1-u1au3μv1-qvIv-u2γh+μh1-u3bu2γh+μhIh(16)-u3bIhIv+u2-u1aShIv.Choosing u1=u2=b/(γh+μh) and u3=1, (16) becomes(17)V′=-u1μhSh0-Sh2Sh-μv1-qv1-R0eIv-bIhIv. Thus, V′ is negative for R0e⩽1. Note also that V′=0 if and only if Sh=Sh0 and Ih=Iv=0. Therefore, the largest invariant set for (1) is the singleton {E10}. Hence, by LaSalle’s invariance principle [34], E10 is globally asymptotically stable when R0e⩽1 and Theorem 3 is valid.Remark 4.
Clearly, the endemic equilibriumE1∗ exists and is unique for R0e>1. This excludes the possibility of occurrence of backward bifurcation. This result is of great epidemiological significance in guiding efforts for disease control as it indicates that R0e=1 is the critical epidemic threshold.To establish the local stability of the endemic equilibrium we evaluate the Jacobian of the system atE1∗, which gives(18)J1∗=-μh+aIv∗0-aSh∗aIv∗-caSh∗0b1-Iv∗-bIh∗-V.The characteristic polynomial of matrix (18) is then given by(19)P2λ=λ3+b2λ2+b1λ+b0,where (20)b2=V2cR0eμh+a+Vc2R0eR0eμh2+a+VaR0eμhvR0e+b+VbR0e2μh+bμh2R0ea+c+μh+abcμh,b1=V2cR0e3μha+c+μh+VbR0e3μh2a+μh+Vc2R0e3μha+μh+VbcR0e2μh2R0e-1+abcR0eμh2R0e-Vμh+bcR0e2μh3,b0=V2c2R0e3μh+VbcR0e2μh2R0e-1+V2ac2R0e2R0e-1+VabcR0eμhR0eR0e-1-1.The coefficient b2>0 and both b1,b0 are nonnegative if and only if R0e>1. Hence, all Routh criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Therefore, the following results holds.Theorem 5.
The endemic equilibriumE1∗=(Sh∗,Ih∗,Iv∗) exists and it is locally asymptotically stable if R0e>1.A global stability result for the endemic equilibriumE1∗ of system (1) is given below.Theorem 6.
IfR0e>1, the endemic equilibrium E1∗=(Sh∗,Ih∗,Iv∗) is globally asymptotically stable.Proof.
LetD1=v1Sh-Sh∗-Sh⁡∗ln⁡(Sh/Sh∗), D2=v2Ih-Ih∗-Ih∗ln⁡(Ih/Ih∗), and D3=v3Iv-Iv∗-Iv∗ln⁡(Iv/Iv∗) be the components of the Lyapunov function(21)USh,Ih,Iv=D1+D2+D3,where v1,v2,v3 are some positive parameters to be chosen later. Differentiating U along the solutions of system (1), we obtain (22)D1′=v1Sh-Sh∗Shμh-μhSh-aShIv=-v1μhSh-Sh∗2Sh+v1Sh-Sh∗ShaSh∗Iv∗-aShIv=-v1μhSh-Sh∗2Sh+Iv∗-ShSh∗Iv-Sh∗ShIv∗-Iv, for v1=1/aSh∗ and μh=μhSh∗-aSh∗Iv∗ at equilibrium. (23)D2′=v2Ih-Ih∗IhaShIv-γh+μhIh=v2Ih-Ih∗IhaShIv-aSh∗Iv∗Ih∗Ih=ShSh∗Iv-IhIh∗Iv∗-Ih∗ShIhSh∗Iv+Iv∗, for v2=v1 and γh+μh=aSh∗Iv∗/Ih∗ at equilibrium. (24)D3′=v3Iv-Iv∗IvbIh1-Iv-μv1-qvIv=bv3Iv-Iv∗IvbIh1-Iv-Ih∗1-Iv∗Iv∗Iv=bv3Iv-Iv∗IvIh-IhIv-IvIh∗1-Iv∗Iv∗,⩽Iv-Iv∗IvIh-IhIv-IvIh∗, for v3=1/b and μv(1-qv)=Ih∗(1-Iv∗)/Iv∗ at equilibrium. Now all together, (25)U′⩽-v1μhSh-Sh∗2Sh-Iv∗IvIv∗+IhIh∗-2-Sh∗ShIv∗-Ih∗ShIhSh∗Iv-Iv∗IvIv∗+IhIh∗-2-2Iv∗+2Iv+Ih-IhIv-IvIh∗-Iv∗IvIh+Iv∗+Iv∗Ih∗,⩽-v1μhSh-Sh∗2Sh-2Iv∗IvIv∗+IhIh∗-2+2Iv1-Iv∗Iv+Ih1-Iv∗Iv+Iv∗Ih∗1-IvIv∗+Iv∗1-Sh∗Sh1+Ih∗Sh2IvIhSh∗2Iv∗. Using the inequality 1-x+ln⁡x⩽0 for x>0 with equality holding if and only if x=1 and the fact that the arithmetic mean is greater than or equal to the geometric mean, we obtain U′(Sh,Ih,Iv)⩽0 for all Sh,Ih,Iv>0. Furthermore, we obtain that U′(Sh,Ih,Iv)=0 holds only when Sh=Sh∗,Ih=Ih∗,Iv=Iv∗ and that E1∗ is the only equilibrium state of these systems on this plane (line). Therefore, by LaSalle’s invariance principle [34], the positive equilibrium E1∗ is globally asymptotically stable.


#### 3.1.2. Epidemic Thresholds, Vertical Infection, and Basic Properties

One of the most important critical thresholds in epidemic models is the basic reproductive number,R0, which is usually found using the next-generation method, as the dominant eigenvalue of thenext-generation matrix [35, 36]. Following the method in [36] we write system (1) consisting only of infectious compartments as the difference between new infection and transfer rates and the resulting Jacobian matrices evaluated at the disease-free equilibrium E0=(1,0,0). For system (1) we have two infected classes, namely, Iv and Ih. It follows that the transmission and the transfer matrices F and F, respectively, are defined as(26)F=qvμvpvhαphvαm0,V=μv00γh+μh. Unlike in host-vector models without vertical transmission, the diagonal elements of the transmission matrix F are nonzero. This stems from the fact that in the presence of vertical transmission there is vector to vector transmission, which completely changes the nature of the basic reproductive number. Thus, the next-generation matrix, K, is then given by(27)K=FV-1=qvpvhαγh+μhphvαmμv0, and the resulting dominant eigenvalue of the spectral radius FV-1, which is the basic reproductive number, R0, is given by(28)R0=12qv+12qv2+4R0H,withR0H=pvhphvα2mμvγh+μh.When there is no vertical transmission, qv=0, as in the case of malaria, R0=R0H is simply the geometric mean of the product of the number of new infections in hosts from one infected vector and the number of new infections in vectors from one infected host, in the limiting case that both populations are fully susceptible. The interpretation and epidemiological significance of R0H are well established. It is easy to see that the transmission of infection is increased with efficiency of vector biting α and probabilities of successful infection transmission, but it is hindered by high mosquito death rates and faster host recovery. The biting rate α appears as α2 because it enters twice in the transmission cycle [31]. An important parameter is the ratio of female mosquitoes to hosts m=Nv/Nh, which is central for disease spread according to model settings. R0 increases with the number (or density) of mosquitoes but decreases with the number (or density) of host population. This results from the asymmetry in the dependence of the vector’s biting rate on the sizes of the host and vector populations such that when there are many more hosts compared to mosquitoes, sustained transmission may be impossible. Therefore, in the absence of vertical transmission, for the infection to successfully spread and invade, the ratio of mosquitoes to hosts needs to be sufficiently large so that double bites are common [37]:(29)NvNh>μvγh+μhpvhphvα2, where the critical ratio is given by mc=μv(γh+μh)/pvhphvα2. Note that each mosquito could infect less than one host on average, and yet R0 could still be more than unity. To elucidate this fact we write R0H as a product of each host type single-step reproductive number; that is,(30)R0H=Rhv×Rvh=phvαmμv×pvhαγh+μh, where Rvh represents the number of new infections in mosquitoes from a single infected host while Rhv represents the number of new infections in hosts from a single infected mosquito. Clearly, R0H can be greater than unity even when one of these reproduction numbers is less than unity, and it can also be less than unity even if one of its components is greater than unity. In Figure 1(a) we depict contours plots corresponding to the overall R0 in (31) along the plane (Rhv,Rvh) in the presence of vertical transmission. Its effects in this asymmetric relationship between R0 and its components is not very pronounced and the ratio of mosquitoes to hosts remains one of the leading factors when there is a large disparity between the sizes of the host and vector populations.Figure 1
Relationships among epidemic thresholds with respect to vertical transmission and ratio of mosquitoes to hosts. (a) The relationship between number of new infections in hosts and number of new infections in vectors. (b) Comparison between basic reproductive numbersR0 and R0e derived from two different methods with respect to vertical infection. (c) Comparison of behaviours among the complete transmission cycle R02 and host type and vector type reproductive numbers with respect to vertical transmission efficiency away from the common epidemic criterion. (d) How the ratio of mosquitoes to hosts m and vertical infection qv contribute to patterns of R0e. Dashed lines indicated contours lines for R0e<1 while continuous lines are for R0e⩾1.
(a) (b) (c) (d)Clearly, the geometric mean is less than the average expected number of new infections per generation. This is the case where the next-generation method fails to produce the correctR0 if transmission between hosts is intermediated by another host (for more discussion about the failure of this method see [21]). Instead, it gives the weighted average lying between the number of new infections each individual produces in the next infection event. Note that if the number of hosts is increased, this deficit is compounded. This can have serious implications for guiding disease control efforts as it fails to provide the actual severity of the infection. Therefore, in this study for application purpose we propose the use of another epidemic threshold. This is R0 from (10) derived in Section 3.1.1, which is hereby referred to as the ‘effective’ reproductive number:(31)R0e=11-qvpvhphvα2mμvγh+μh,as (1) it satisfies the property that the endemic equilibrium E1∗ only persists if R0e is greater than unity and (2) the endemic equilibrium exists without occurrence of backward bifurcation, meaning that R0e=1 is the correct critical value. At this point, it is of particular interest to establish the critical ratio of mosquitoes to hosts in the settings of vertical transmission. From (31) we obtain(32)NvNh>1-qvμvγh+μhpvhphvα2, that is, a new critical ratio mc∗=(1-qv)mc<mc for 0<qv<1. This result indicates that, in the presence of vertical transmission, the critical ratio of female mosquitoes to hosts reduces at a rate proportional to vertical transmission efficiency. This indicates an earlier occurrence of the point where an epidemic is just possible as above this level the equilibrium prevalence is expected to rapidly increase to its asymptotic value.Others are epidemic thresholds that provide a direct measure of the control effort required for disease eradication. These measures are known as “type” reproductive numbers and are determined based on the next-generation matrix [38, 39] such that for our case we have host and vector type reproductive numbers T1h and T1v, respectively. Note that in the absence of vertical transmission the host and vector type reproductive numbers are both equal to the square of the basic reproductive number, R0 in (41), which turns out to be R0e in (31) for qv=0. In the presence of vertical infection, R02 gives the expected number of secondary infections after one average, complete (host-vector-host or vector-host-vector) transmission cycle but does not correspond to a specific population type [23]. For the case of host population, one infected host leads to some secondary host infections in the next host-vector-host transmission cycle. This results exactly from the horizontal transmission mechanism and it is given by R0H. However, further secondary host infections may also occur after any number of vector-vector transmission cycles as a result of transovarial transmission. From the expression of R0 in (31) it can be seen that new infections resulting from vector-vector transmission cycle are given by 1/(1-qv); hence the host type reproductive number is then given by(33)T1h=11-qvR0H=R0e.From (33), R0H=R0e if 1/(1-qv)=1; that is, qv∗=0 which is the boundary condition. An important question would be, what is the proportion of vertical infection that doubles R0H? For 1/(1-qv)=2 we obtain qv∗=0.5, which means that if about 50% of infected mosquitoes produce infected offspring, we expect T1h to double. Note that all secondary infections resulting from a single infected vector must occur in the next vector-host-vector transmission cycle; hence, the vector type reproductive number takes the form(34)T1v=R0H+qv.For details on the derivation of both types of reproductive numbers, see Supplementary Materials, Section A-1. A number of methods for deriving R0 exist and in [21] it is demonstrated that the resulting R0 is not unique but they may have a common threshold at R0=1. Hence, their behaviour below and above unity may show significant differences. In Figure 1© we examine how vertical transmission efficiency contributes to the qualitative behaviour of R02,T1h, and T1v above the unity threshold criterion. An important characteristic pattern of their relationship is that below a certain percentage of vertical transmission efficiency the three reproductive numbers are undistinguishable but beyond the relationship between vertical infection and host type reproductive number becomes exponential. This behaviour results from the fact that in each generation, the number of host infections is proportional to the number of infected mosquitoes, which decreases proportionally to the vertical infection efficiency. Hence, for low vertical infection efficiency, possibly the first generations make a significant contribution to subsequent transmission cycles (for more discussion see [23, 27]).In Figure1(b) we compare the basic reproductive number R0 derived through the next-generation method with the one derived from the endemic equilibrium R0e with respect to vertical transmission efficiency. Clearly, the two give different relationships and it is easy to see that the effect of vertical infection is completely diluted on R0. Hence, this measure if used for disease control effort is more likely to provide misleading guidance. The relationship of the basic reproductive number R0e with mosquito to host ratio m and vertical infection qv shows some interesting patterns. For m>1 the relationship between R0e and qv is exponential as described in Figure 1(c); however for m<1 it is linear rendering qv having little effect if not negligible for qv∈(0,0.8) (see Figure 1(d) dashed lines). This suggests that when there are more hosts compared to mosquitoes, vertical infection efficiency becomes completely negligible as in the model settings, and if m<1, sustained host-vector transmissions may be impossible. On the other hand, we observe that for qv⋙0.8 and even for m<1, R0e becomes much larger highlighting the importance of the interplay between vertical infection and the ratio of mosquitoes to hosts.Remark 7.
In the settings of vector-borne diseases both host type reproductive numberT1h or (R0e) and vector type reproductive number T1v are critical for providing useful informative indicator for guiding disease control efforts.


### 3.2. Effects of Ratio of Mosquitoes to Hosts and Vertical Transmission on Disease Prevalence

To fully capture effects of vertical transmission efficiency on vector-borne transmission models, we explore how this parameter influences the disease spread patterns.Results are shown in Figure2. The left figure explores how the first peak of an epidemic denoted by max(Ih(t)) varies with vertical transmission efficiency, qv. The results show that max(Ih(t)) increases linearly with qv. The right-hand figure plots the changes on time of occurrence of each epidemic peak at each value of vertical transmission efficiency. The results suggest that as vertical transmission efficiency increases the time position of each epidemic peak decreases linearly. This suggests that if vertical infection efficiency is increased, the first peak of an epidemic is likely to be attained earlier than the expected time. This suggests that while vertical transmission efficiency increases the epidemic size of an outbreak, it linearly decreases the time of the peak of an outbreak. The lower figures, (c) and (d), describes how the prevalence of infected hosts, Ih, varies with the ratio of female mosquitoes to hosts for different values of vertical transmission efficiency. Figure 2(c) indicates that the prevalence of the disease saturates at larger values of the ratio at earlier stage of the initial epidemic while Figure 2(d) indicates that the prevalence of the disease saturates at lower values of the ratio at later stage of the initial epidemic. This result stems from the fact that at earlier stage of the initial epidemic the size of the vector population is still very high compared to later stages of evolution of the epidemic.Figure 2
The left figure (a) explores effects of vertical transmission on the prevalence of the disease on hosts, where (Ih(qv)) is plotted at max(Ih(t)) for every qv. The right figure (b) explores how the first peak of an epidemic (max(Ih(t))) varies with vertical transmission efficiency (qv). The lower figures represent the dynamics of Ih(t) at time position (78) corresponding to day 27 (figure c) and at time position (150) corresponding to day 50 (figure d) for different values of vertical infection.
(a) (b) (c) (d)


### 3.3. Analysis of the Model with Vector Stages and Vertical Transmission

#### 3.3.1. Model Equilibria

Model equilibria are obtained by setting the right-hand side of system (2) equal to zero. The model has three equilibria, two disease-free and one endemic equilibrium. Details about the derivation of the components of the endemic equilibrium are given in Section A-2 of Supplementary Materials. The first trivial equilibrium is both mosquito-free and disease-free:(35)X0=Ls0,Li0,Sv0,Iv0,Sh0,Ih0=0,0,0,0,1,0, which corresponds only to the presence of hosts (this can be livestock for the case of RVF). However, this equilibrium is not biologically feasible or relevant in particular in tropical regions where vector-borne diseases such as Dengue and RVF are endemic. The second equilibrium corresponds to the coexistence of both vectors and hosts but without infection. This is the disease-free equilibrium:(36)X1=Ls1,Li1,Sv1,Iv1,Sh1,Ih1=l∗,0,P∗,0,1,0, which can be used for computation of the basic reproduction number when using the next-generation method [36]. The nonzero components are(37)l∗=δr1-Q0-1δr+θ+μL,P∗=θ1-Q0-1θ+δμv,withQ0=rθμvθ+μL. Hence, X1 is biologically significant if Q0>1; that is, the mosquito population exists and can establish itself if on average each adult female mosquito produces at least one female mosquito at oviposition rate r during her life time period 1/μv after successfully surviving the aquatic stage at probability θ/(θ+μL).The endemic equilibrium which represents the prevalence of the disease in both host and vector populations is given by(38)X∗=Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗, where(39)Ls∗=l∗-Li∗,Sv∗=P∗-Iv∗,Iv∗=μh1-Sh∗λhvSh∗,Sh∗=1-γh+μhμhIh∗,(40)Ih∗=1-qvμvγh+μhR0c-1λvhλhvP∗γh+μh/μh+λvh/1-qvμv,where(41)R0c=λvhλhvP∗1-qvμvγh+μh.Similar to the discussion in Section 3.1.2 the above epidemic threshold in (41) is the “effective” or “control” reproductive number since it satisfies the property that the endemic equilibrium E∗ only persists if R0c is greater than unity and there is no backward bifurcation.Remark 8.
The relationship between the endemic equilibrium andR0 which depends on the vector population threshold indicates that the existence of this equilibrium is governed by the norms of the vector population threshold. Therefore, this vector population threshold is an important parameter in vector control efforts as it provides a means for identifying key factors for reducing the vector population.In the following section we deriveR0 using the next-generation method.


#### 3.3.2. The Basic Reproductive Number and Other Control Thresholds

The basic reproductive number,R0, is a concept of great epidemiological significance. Following the method of [36], we write system (2) considering only disease compartments, Li,Iv and Ih as x˙=F-V with xT=(Li,Iv,Ih), where(42)F=δr1-LqvIvθδ1-NvLi+λvhIhSvλhvIvSh,V=θ+μLLiμvIvγh+μhIh. Evaluating the Jacobian of above matrices at the disease-free equilibrium X1=(l∗,0,P∗,0,1,0), the basic reproduction number is then the spectral radius of the next-generation matrix; that is, R0=ρ(F(X1)V(X1)-1),(43)R0=12qv+12qv2+4R0H,withR0H=λvhλhvP∗μvγh+μh.Following the discussion in Section 3.1.2 this epidemic threshold (43) can lead to an underestimate of the control effort required if used to guide disease control strategies. To overcome this inconsistency we use type reproductive numbers to estimate efforts required to provide informative indicators in control of vector-borne diseases. To avoid repetition we omit the steps required in the derivation process, such that the host type reproductive number is(44)T1h=R0H1-qv=R0c,and the vector type reproductive number is then given by(45)T1v=R0H+qv.In the majority of mosquito-borne diseases, treatment or adequate vaccine is not likely to be feasible, either due to financial constraints or nonavailability, in particular in African communities with low resilience to economic challenges. In such cases the main preventive measures are individual protection against mosquito bites and the control of the proliferation of both the larval and adult mosquitoes. This can be achieved by lowering appropriate epidemic thresholds below unity, that is, forcing the reproductive number of the pathogen below 1. In the context of our model this implies reducing either of or both type reproductive numbers (T1h and/or T1v) below unity. Hence, we define the eradication effort to be the percentage reduction in the mosquito population size required to prevent endemic transmission of the pathogen [23], such that if using T1h the eradication effort is 100(1-1/T1h)=100(R0H+qv-1)/R0H and for T1v the effort required to eradicate the disease is 100(1-1/T1v)=100(R0H+qv-1)/(R0H+qv). Therefore, vertical infection at qv% can be responsible for at most qv% of the required eradication effort, provided R0H>1.


#### 3.3.3. Stability of the Host-Only System

This fixed pointX0=(Ls0,Li0,Sv0,Iv0,Sh0,Ih0)=(0,0,0,0,1,0) is an equilibrium point of the host-only system given by (46)dShdt=μh1-Sh, which is local and globally asymptotically stable for all t>0 as the derivative μh1-Sh′ is negative.


#### 3.3.4. Local Stability of the Disease-Free Equilibrium

The local stability of the disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is established by analysing the eigenvalues of the Jacobian matrix at X1:(47)JE1=-δrP∗-θ+μLδrP∗δr1-l∗δr1-l∗1-qv000-θ+μL0δr1-l∗qv00θδ1-P∗0-θδl∗-μv-θδl∗0-λvhP∗0θδ1-P∗0-μv0λvhP∗000-λhv-μh0000λhv0-γh+μh.Direct computation shows that λ1=-μh is an eigenvalue of matrix (47) and the remaining are solutions of equation(48)λ5+c4λ4+c3λ3+c2λ2+c1λ+c0=0where(49)c4=δrP∗+2θ+μL+γh+μh+θδl∗+2μv>0,c3=1-qvγh+μhμv1-R0c+δrP∗γh+μh+θl∗δ+2μv+θ+μL+μv2+θδl∗2θ+μL+γh+μh+μv+2γh+μhθ+μL+μv+1-qvμvθ+μL-γh+μh>0,c2=1-qvγh+μhμv1-R0cδrP∗+2θ+μL+θl∗δ+μv+1-qvμvθ+μL-γh+μhδrP∗+θ+μL+θl∗δ+μv+γh+μhθ+μL+μv2+δrP∗θ+μLγh+μh+θl∗δ+μv+γh+μhθl∗δ+2μv+μvθl∗δ+μv+θδl∗θ+μLθ+μL+2γh+μh+μv+θδl∗γh+μhμv>0,c1=1-qvγh+μhμv1-R0cδrP∗θ+μL+θl∗δ+μv+θ+μL2+θ+μLμv+2θ+μLθl∗δ+δrP∗θl∗δ+μv+θ+μLθl∗δ1-qvμvθ+μL-γh+μh+γh+μhμv+θ+μLγh+μhθl∗δδrP∗+θ+μL>0,c0=1-qvθ+μLγh+μhμv1-R0cδrP∗θl∗δ+μv+θl∗δμv.By Descartes’ rule of signs, since c4>0, c3>0, c2>0, and c1>0, then all roots of (48) are negative or have a negative real part for c0>0. The coefficient c0 is nonnegative if and only if R0c<1; hence the following result holds.Theorem 9.
The disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is locally asymptotically stable whenever R0c<1.


#### 3.3.5. Global Stability of the Disease-Free Equilibrium

To establish the global asymptotic stability (GAS) ofX1 we use results obtained by Kamgang and Sallet [40], which are the extension of some results in [36]. Around the disease-free equilibrium system (2) can be written as(50)x˙S=A1xxS-xE1,S+A12xxIx˙I=A2xxI where xS is the vector representing disease-free compartments (Ls,Sv,Sh) and the vector xI represents the state of infected compartments (Li,Iv,Ih). This requires rewriting system equation around X1, as follows: (51)dShdt=μh-λhvIvSh-μhSh=-μhSh-Sh1-λhvIvSh,dLsdt=δr1-LsSv+1-qvIv-δrLiSv+1-qvIv-θ+μLLs,=-δrSv+θ+μLLs-Ls1+δr1-Ls1Sv-Sv1-δrSvLi+δr1-qv1-LIv.Note that the latter arises from (θ+μL)Ls1=δr(1-Ls1)Sv1. Then, the following matrices are obtained,(52)A1x=-δrSv+θ+μLδr1-Ls10θδ1-Nv-μv000-μh,A12x=-δrSvδr1-qv1-L000-λvhSv0-λvhSv0,(53)A2x=-θ+μLδr1-Lqv0θδ1-Nv-μvλvhSv0λhv-γh+μh. From basic matrix properties and direct algebraic computation it follows that all eigenvalues of matrix A1 are real and negative and both A1 and A2 are Metzler matrices. Hence, system x˙S=A1(x)(xS-xX1,S) is GAS at the disease-free equilibrium xX1,S. To establish the global stability of overall system (50) at X1 conditions of the following theorem in [40] must be satisfied.Theorem 10.
LetΦ⊂U=R+3×R+3. System (50) is of class C1, defined on U if(1)
U is positively invariant relative to (50).(2)
The systemx˙S=A1(x)(xS-xDFE,S) is GAS at xDFE,S.(3)
For anyx∈Φ, matrix A2(x) is Metzler irreducible.(4)
There exists a matrixA¯2, which is an upper bound of the set M=A2(x)∈M3(R)∣x∈Φ¯, with the property that if A¯2∈M, for any x¯∈Φ¯, such that A2(x¯)=A¯2, then x¯∈R3×{0}.(5)
The stability modulus ofA¯2,α(A¯2)=maxλ∈Sp(A2)Re(λ), satisfies α(A¯2)≤0.
Then,X1 is GAS in Φ¯.Proof.
Clearly, conditions(1)-(3) of the theorem have been satisfied. For all x∈Φ, A2(x) is irreducible because (I+A2x)2>0. An upper bound of the set of matrices M, which is the matrix A¯2, is given by matrix A2(x¯), where x¯=(L¯s1,0,S¯v1,0,1,0)∈R3×{0}, with L¯s1=K1 and S¯v1=K2. Similarly matrix A¯2 is irreducible. Recall that from the Perron-Frobenius theorem for an irreducible matrix you get that one of the matrix eigenvalues is positive and greater than or equal to all others, that is, the dominant eigenvalue. Matrix A2 is exactly the matrix used to compute the basic reproductive number, i.e., the dominant eigenvalue. For more details or proof in general settings see [40].Now conditions(1)-(4) have been verified. To check the last condition, we make use of the following Lemma [40].Lemma 11.
LetH be a square Metzler matrix written in block form H=(ABCD), with A and D square matrices. H is Metzler stable if and only if matrices A and D-CA-1B are Metzler stable.MatrixA2(x) in block matrices takes the following components:(54)A=-θ+μL,B=δr1-Lqv0,C=θδ1-Nv0T,D=-μvλvhP∗λhv-γh+μh. Clearly, A is a Metzler stable matrix and(55)D-CA-1B=-μv1-qvλvhP∗λhv-γh+μh is also Metzler stable if μv(1-qv)(γh+μh)-λhvλvhP∗⩾0, that is, 1-R0c⩾0. Therefore, from Theorem 10 and Lemma 11 the following result holds.Theorem 12.
ForQ0>1 there exists a unique disease-free equilibrium X1=(l∗,0,P∗,0,1,0), which is globally asymptotically stable whenever R0c⩽1.Remark 13.
The above result is of great epidemiological importance. It highlights two fundamental indicators:(1) Although vector control is central for control of vector-borne diseases, it does not mean eliminating all the vectors. Note that the effect of this include both infected and noninfected mosquitoes. (2) Disease eradication efforts are independent of the initial sizes of both vector and host populations; however, the ratio between the two populations is a key factor.


#### 3.3.6. Local Stability of the Endemic Equilibrium

Results of the local stability of the disease-free equilibriumX1 suggest that for R0c=1 the Jacobian matrix (47) has zero eigenvalue while the remaining eigenvalues are negative or have a negative real part. Since the algebraic computation involved when establishing the stability of the endemic equilibrium through linearization is quite extensive, we employ the centre manifold theory [41]. This theory is used to examine existence of backward or forward bifurcation. The bifurcation occurs at R0c=1, and choosing qv=qv∗ as a bifurcation parameter, then(56)qv∗=1-λhvλvhP∗μvγh+μh,⇒R0H⩽1.The Jacobian matrix J(X1,qv∗) is the same as matrix (47); hence, J(X1,qv∗) has a simple zero eigenvalue λ=0 when R0c=1 and others given as λ=-μh and roots of the polynomial λ4+c4λ3+c3λ2+c2λ+c1=0, where c1,c2,c3,c4 are as defined in Section 3.3.4. Let ω=(ω1,ω2,ω3,ω4,ω5,ω6) be the right eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are derived by solving J(X1,qv∗)×ω=0, which gives(57)ω=F∗ω4,H∗ω4,G∗ω4,ω4,-λhvμhω4,λhvγh+μh where (58)F∗=rθμhQ0-1-qv∗+1-qv∗μhγh+μhθ+μLδrθ+δμvrθ2μh1-Q0-1δr+θ+μL,G∗=1-qv∗μvθθ+μLγh+μhδr+θ+μL-μhθ+μLQ0v-1+μhμvθ+μL2θ+δμv1-Q0qv∗rθθ+δμv1-Q0-1,H∗=δrqv∗θ+δμvθδr+θ+μL,Q0v=rθ1-qv∗μvθ+μL.Clearly, F∗ and G∗ are nonnegative if and only if qv∗⩽1/Q0, Q0v>1, and Q0>1. Note that Q0v gives the average number of emerged infected female mosquitoes produced by one infected female mosquito. Hence, for Q0v>1 infected mosquitoes will persist in the population, leading to disease persistence. Let us denote υ=(υ1,υ2,υ3,υ4,υ5,υ6) as the left eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are obtained when solving [J(X1,qv∗)]T×ω=0, (T represents transpose), which gives(59)υ=0,M∗υ4,0,υ4,0,λvhP∗γh+μhυ4 where M∗=θμv(δr+θ+μL)/δr(θ+μL)(θ+δμv).Computation of Bifurcation Parameters a ∗ and b ∗. Because υ1=υ3=υ5=0 the associated nonzero partial derivatives of f1,f3,f5 are not considered and the remaining nonzero partial derivatives of system f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) are given by ∂2f2/∂x1∂x4=∂2f2/∂x2∂x4=-δrqv∗, ∂2f4/∂x2∂x4=∂2f4/∂x3∂x2=-θ/δ, ∂2f4/∂x3∂x6=λvh, and ∂2f6/∂x4∂x5=λhv.Thus,a∗=2υ2ω1ω4(∂2f2/∂x1∂x4)+2υ2ω2ω4(∂2f2/∂x2∂x4)+2υ4ω2ω4(∂2f4/∂x2∂x4)+2υ4ω3ω2(∂2f4/∂x3∂x2)+2υ4ω3ω6(∂2f4/∂x3∂x6)+2υ6ω4ω5(∂2f6/∂x4∂x5), which after some algebraic simplification gives(60)a∗=-2υ4ω42δrqv∗F∗M∗+δrqv∗H∗M∗+θδH∗+G∗H∗+1-qv∗μvλvhμh-G∗P∗.Given that the vital dynamics of the host are much more slower λvh/μh⋙G∗/P∗ and using the property that υω=1, then a∗<0 for ω4>0. For b∗ the nonzero partial derivative of f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) is ∂2f2/∂x4∂qv∗=δr(1-l∗), such that b∗=υ2ω4(∂2f2/∂x4∂qv∗)=μvυ2ω4. Using the property that υω=1 we obtain that b∗>0. Since a∗<0 for ω4>0 and b∗>0, the model (2) exhibits a forward bifurcation at R0c=1. Therefore, the following result holds.Theorem 14.
The endemic equilibrium pointX∗=(Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗) exists and is locally asymptotically stable for R0c>1.Remark 15.
The conditionQ0v>1 is an indication that vertical transmission is vital for disease long-term persistence. A particular example is the case of Dengue fever and RVF. Upon failure of this condition, it is much more likely that results of Theorem 14 do not hold. Another important feature of Q0v is that vertical infection is proportional to Q0, and hence vector control is still a viable control strategy even for vector-borne disease with transovarial transmission.The above stability results and those in Theorems9, 10, 12, and 14 can be summarized in Table 1 and in the bifurcation diagram in Figure 3.Table 1
Existence and stability of model equilibria.
X 0 X 1 if (0<qv<1) X ∗ if (0<qv<1) I. Q 0 < 1   &   R 0 c < 1 GS DNE DNE II. Q 0 > 1   &   R 0 c < 1 DNE GS DNE III. Q 0 > 1   &   R 0 c > 1 DNE US LS IV. Q 0 < 1   &   R 0 c > 1 US DNE DNE DNE: does not exist, US: unstable, LS: locally stable, GS: globally stable.Figure 3
The bifurcation diagram in (Q0,R0c) plane. The stability of equilibria in each region is given in Table 1.


### 3.4. Numerical Results of the Model with Vector Stages and Vertical Transmission

Numerical analysis using reasonable parameter values for RVF is carried out. In Figure4 we plot time series of infected larval, adult mosquitoes and hosts, and their respective phase portraits. The latter are plotted for t∈[050] years for observing long-term dynamics. The parameter values used give R0c>1 and existence of a stable endemic equilibrium. These results not only confirm findings by analytical analysis (Theorems 9 and 12) but suggest that this equilibrium is also global asymptotic stable for the chosen parameter values. This also implies that the endemic equilibrium is stable not only close to the bifurcation point but also for R0c beyond unity (see inner figures, Figures 4(a)–4(c)).Figure 4
Time series and phase diagrams plot. The inner figures describe the temporal evolution ofLi,Iv,Ih after the transient period for R0c=2.1002, while for the outer figures R0c=1.0012. Parameter values used for R0c=1.0012 are μv=1/25,μh=1/(8∗360),phv=0.14,pvh=0.35,α=256/360, γh=1/4,m=1.5,qv=0.2,θ=1/10,μL=1/18,r=1/12,δ=0.4 and parameter values used for R0c=2.1002 are given in Table 1 of Section A-3 of Supplementary Materials.
(a) (b) (c) (d) (e) (f)


## 3.1. Analysis of the Basic Host-Vector Model

### 3.1.1. Model Equilibria and Stability Analysis

The basic host-vector model with vertical transmission exhibits two equilibria, namely, the disease-freeE10 and the endemic E1∗, respectively. At the disease-free equilibrium,(8)E10=Sh0,Ih0,Iv0=1,0,0,both vector and host populations persist but with no disease. The prevalence of the disease is denoted by(9)E1∗=Sh∗,Ih∗,Iv∗=a+μhR0ea+μhR0h,μhR0e-1cR0e+bμh/V,μhR0e-1a+μhR0e,where a=phvαm, b=pvhα, c=γh+μh, V=μv(1-qv), and(10)R0e=abVc=phvpvhα2mμv1-qvγh+μh.The Jacobian matrix of system (1) at E10 is given by(11)J1′=-μh0-a0-ca0b-V,such that the characteristic polynomial of matrix (11) is then given as(12)P1λ=λ3+a2λ2+a1λ+a0,with a2=V+c+μh, a1=Vc(1-R0e)+(V+c)μh, and a0=μhVc(1-R0e). The coefficient a2>0 and both a1,a0 are nonnegative if and only if R0e<1. Hence, all Routh stability criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Furthermore, for R0e=1, (9) becomes the disease-free equilibrium. Therefore, the following result holds.Theorem 3.
The disease-free equilibriumE10=(1,0,0) exists and it is globally asymptotically stable if R0e⩽1.Alternatively, the global stability of the disease-free equilibriumE10 can be established using the following Lyapunov function:(13)VSh,Ih,Iv=u1Sh-Sh0-Sh0lnShSh0+u2Ih+u3Ivwhere u1,u2,u3 are some positive constants. Calculating the derivative of V along the solutions of system (1), we obtain(14)V′=u1Sh-Sh0ShμhSh0-Sh-aShIv+u2aShIv-γh+μhIh+u3bIh1-Iv-μv1-qvIv,(15)=-u1μhSh0-Sh2Sh-u3μv1-qv1-u1au3μv1-qvIv-u2γh+μh1-u3bu2γh+μhIh(16)-u3bIhIv+u2-u1aShIv.Choosing u1=u2=b/(γh+μh) and u3=1, (16) becomes(17)V′=-u1μhSh0-Sh2Sh-μv1-qv1-R0eIv-bIhIv. Thus, V′ is negative for R0e⩽1. Note also that V′=0 if and only if Sh=Sh0 and Ih=Iv=0. Therefore, the largest invariant set for (1) is the singleton {E10}. Hence, by LaSalle’s invariance principle [34], E10 is globally asymptotically stable when R0e⩽1 and Theorem 3 is valid.Remark 4.
Clearly, the endemic equilibriumE1∗ exists and is unique for R0e>1. This excludes the possibility of occurrence of backward bifurcation. This result is of great epidemiological significance in guiding efforts for disease control as it indicates that R0e=1 is the critical epidemic threshold.To establish the local stability of the endemic equilibrium we evaluate the Jacobian of the system atE1∗, which gives(18)J1∗=-μh+aIv∗0-aSh∗aIv∗-caSh∗0b1-Iv∗-bIh∗-V.The characteristic polynomial of matrix (18) is then given by(19)P2λ=λ3+b2λ2+b1λ+b0,where (20)b2=V2cR0eμh+a+Vc2R0eR0eμh2+a+VaR0eμhvR0e+b+VbR0e2μh+bμh2R0ea+c+μh+abcμh,b1=V2cR0e3μha+c+μh+VbR0e3μh2a+μh+Vc2R0e3μha+μh+VbcR0e2μh2R0e-1+abcR0eμh2R0e-Vμh+bcR0e2μh3,b0=V2c2R0e3μh+VbcR0e2μh2R0e-1+V2ac2R0e2R0e-1+VabcR0eμhR0eR0e-1-1.The coefficient b2>0 and both b1,b0 are nonnegative if and only if R0e>1. Hence, all Routh criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Therefore, the following results holds.Theorem 5.
The endemic equilibriumE1∗=(Sh∗,Ih∗,Iv∗) exists and it is locally asymptotically stable if R0e>1.A global stability result for the endemic equilibriumE1∗ of system (1) is given below.Theorem 6.
IfR0e>1, the endemic equilibrium E1∗=(Sh∗,Ih∗,Iv∗) is globally asymptotically stable.Proof.
LetD1=v1Sh-Sh∗-Sh⁡∗ln⁡(Sh/Sh∗), D2=v2Ih-Ih∗-Ih∗ln⁡(Ih/Ih∗), and D3=v3Iv-Iv∗-Iv∗ln⁡(Iv/Iv∗) be the components of the Lyapunov function(21)USh,Ih,Iv=D1+D2+D3,where v1,v2,v3 are some positive parameters to be chosen later. Differentiating U along the solutions of system (1), we obtain (22)D1′=v1Sh-Sh∗Shμh-μhSh-aShIv=-v1μhSh-Sh∗2Sh+v1Sh-Sh∗ShaSh∗Iv∗-aShIv=-v1μhSh-Sh∗2Sh+Iv∗-ShSh∗Iv-Sh∗ShIv∗-Iv, for v1=1/aSh∗ and μh=μhSh∗-aSh∗Iv∗ at equilibrium. (23)D2′=v2Ih-Ih∗IhaShIv-γh+μhIh=v2Ih-Ih∗IhaShIv-aSh∗Iv∗Ih∗Ih=ShSh∗Iv-IhIh∗Iv∗-Ih∗ShIhSh∗Iv+Iv∗, for v2=v1 and γh+μh=aSh∗Iv∗/Ih∗ at equilibrium. (24)D3′=v3Iv-Iv∗IvbIh1-Iv-μv1-qvIv=bv3Iv-Iv∗IvbIh1-Iv-Ih∗1-Iv∗Iv∗Iv=bv3Iv-Iv∗IvIh-IhIv-IvIh∗1-Iv∗Iv∗,⩽Iv-Iv∗IvIh-IhIv-IvIh∗, for v3=1/b and μv(1-qv)=Ih∗(1-Iv∗)/Iv∗ at equilibrium. Now all together, (25)U′⩽-v1μhSh-Sh∗2Sh-Iv∗IvIv∗+IhIh∗-2-Sh∗ShIv∗-Ih∗ShIhSh∗Iv-Iv∗IvIv∗+IhIh∗-2-2Iv∗+2Iv+Ih-IhIv-IvIh∗-Iv∗IvIh+Iv∗+Iv∗Ih∗,⩽-v1μhSh-Sh∗2Sh-2Iv∗IvIv∗+IhIh∗-2+2Iv1-Iv∗Iv+Ih1-Iv∗Iv+Iv∗Ih∗1-IvIv∗+Iv∗1-Sh∗Sh1+Ih∗Sh2IvIhSh∗2Iv∗. Using the inequality 1-x+ln⁡x⩽0 for x>0 with equality holding if and only if x=1 and the fact that the arithmetic mean is greater than or equal to the geometric mean, we obtain U′(Sh,Ih,Iv)⩽0 for all Sh,Ih,Iv>0. Furthermore, we obtain that U′(Sh,Ih,Iv)=0 holds only when Sh=Sh∗,Ih=Ih∗,Iv=Iv∗ and that E1∗ is the only equilibrium state of these systems on this plane (line). Therefore, by LaSalle’s invariance principle [34], the positive equilibrium E1∗ is globally asymptotically stable.


### 3.1.2. Epidemic Thresholds, Vertical Infection, and Basic Properties

One of the most important critical thresholds in epidemic models is the basic reproductive number,R0, which is usually found using the next-generation method, as the dominant eigenvalue of thenext-generation matrix [35, 36]. Following the method in [36] we write system (1) consisting only of infectious compartments as the difference between new infection and transfer rates and the resulting Jacobian matrices evaluated at the disease-free equilibrium E0=(1,0,0). For system (1) we have two infected classes, namely, Iv and Ih. It follows that the transmission and the transfer matrices F and F, respectively, are defined as(26)F=qvμvpvhαphvαm0,V=μv00γh+μh. Unlike in host-vector models without vertical transmission, the diagonal elements of the transmission matrix F are nonzero. This stems from the fact that in the presence of vertical transmission there is vector to vector transmission, which completely changes the nature of the basic reproductive number. Thus, the next-generation matrix, K, is then given by(27)K=FV-1=qvpvhαγh+μhphvαmμv0, and the resulting dominant eigenvalue of the spectral radius FV-1, which is the basic reproductive number, R0, is given by(28)R0=12qv+12qv2+4R0H,withR0H=pvhphvα2mμvγh+μh.When there is no vertical transmission, qv=0, as in the case of malaria, R0=R0H is simply the geometric mean of the product of the number of new infections in hosts from one infected vector and the number of new infections in vectors from one infected host, in the limiting case that both populations are fully susceptible. The interpretation and epidemiological significance of R0H are well established. It is easy to see that the transmission of infection is increased with efficiency of vector biting α and probabilities of successful infection transmission, but it is hindered by high mosquito death rates and faster host recovery. The biting rate α appears as α2 because it enters twice in the transmission cycle [31]. An important parameter is the ratio of female mosquitoes to hosts m=Nv/Nh, which is central for disease spread according to model settings. R0 increases with the number (or density) of mosquitoes but decreases with the number (or density) of host population. This results from the asymmetry in the dependence of the vector’s biting rate on the sizes of the host and vector populations such that when there are many more hosts compared to mosquitoes, sustained transmission may be impossible. Therefore, in the absence of vertical transmission, for the infection to successfully spread and invade, the ratio of mosquitoes to hosts needs to be sufficiently large so that double bites are common [37]:(29)NvNh>μvγh+μhpvhphvα2, where the critical ratio is given by mc=μv(γh+μh)/pvhphvα2. Note that each mosquito could infect less than one host on average, and yet R0 could still be more than unity. To elucidate this fact we write R0H as a product of each host type single-step reproductive number; that is,(30)R0H=Rhv×Rvh=phvαmμv×pvhαγh+μh, where Rvh represents the number of new infections in mosquitoes from a single infected host while Rhv represents the number of new infections in hosts from a single infected mosquito. Clearly, R0H can be greater than unity even when one of these reproduction numbers is less than unity, and it can also be less than unity even if one of its components is greater than unity. In Figure 1(a) we depict contours plots corresponding to the overall R0 in (31) along the plane (Rhv,Rvh) in the presence of vertical transmission. Its effects in this asymmetric relationship between R0 and its components is not very pronounced and the ratio of mosquitoes to hosts remains one of the leading factors when there is a large disparity between the sizes of the host and vector populations.Figure 1
Relationships among epidemic thresholds with respect to vertical transmission and ratio of mosquitoes to hosts. (a) The relationship between number of new infections in hosts and number of new infections in vectors. (b) Comparison between basic reproductive numbersR0 and R0e derived from two different methods with respect to vertical infection. (c) Comparison of behaviours among the complete transmission cycle R02 and host type and vector type reproductive numbers with respect to vertical transmission efficiency away from the common epidemic criterion. (d) How the ratio of mosquitoes to hosts m and vertical infection qv contribute to patterns of R0e. Dashed lines indicated contours lines for R0e<1 while continuous lines are for R0e⩾1.
(a) (b) (c) (d)Clearly, the geometric mean is less than the average expected number of new infections per generation. This is the case where the next-generation method fails to produce the correctR0 if transmission between hosts is intermediated by another host (for more discussion about the failure of this method see [21]). Instead, it gives the weighted average lying between the number of new infections each individual produces in the next infection event. Note that if the number of hosts is increased, this deficit is compounded. This can have serious implications for guiding disease control efforts as it fails to provide the actual severity of the infection. Therefore, in this study for application purpose we propose the use of another epidemic threshold. This is R0 from (10) derived in Section 3.1.1, which is hereby referred to as the ‘effective’ reproductive number:(31)R0e=11-qvpvhphvα2mμvγh+μh,as (1) it satisfies the property that the endemic equilibrium E1∗ only persists if R0e is greater than unity and (2) the endemic equilibrium exists without occurrence of backward bifurcation, meaning that R0e=1 is the correct critical value. At this point, it is of particular interest to establish the critical ratio of mosquitoes to hosts in the settings of vertical transmission. From (31) we obtain(32)NvNh>1-qvμvγh+μhpvhphvα2, that is, a new critical ratio mc∗=(1-qv)mc<mc for 0<qv<1. This result indicates that, in the presence of vertical transmission, the critical ratio of female mosquitoes to hosts reduces at a rate proportional to vertical transmission efficiency. This indicates an earlier occurrence of the point where an epidemic is just possible as above this level the equilibrium prevalence is expected to rapidly increase to its asymptotic value.Others are epidemic thresholds that provide a direct measure of the control effort required for disease eradication. These measures are known as “type” reproductive numbers and are determined based on the next-generation matrix [38, 39] such that for our case we have host and vector type reproductive numbers T1h and T1v, respectively. Note that in the absence of vertical transmission the host and vector type reproductive numbers are both equal to the square of the basic reproductive number, R0 in (41), which turns out to be R0e in (31) for qv=0. In the presence of vertical infection, R02 gives the expected number of secondary infections after one average, complete (host-vector-host or vector-host-vector) transmission cycle but does not correspond to a specific population type [23]. For the case of host population, one infected host leads to some secondary host infections in the next host-vector-host transmission cycle. This results exactly from the horizontal transmission mechanism and it is given by R0H. However, further secondary host infections may also occur after any number of vector-vector transmission cycles as a result of transovarial transmission. From the expression of R0 in (31) it can be seen that new infections resulting from vector-vector transmission cycle are given by 1/(1-qv); hence the host type reproductive number is then given by(33)T1h=11-qvR0H=R0e.From (33), R0H=R0e if 1/(1-qv)=1; that is, qv∗=0 which is the boundary condition. An important question would be, what is the proportion of vertical infection that doubles R0H? For 1/(1-qv)=2 we obtain qv∗=0.5, which means that if about 50% of infected mosquitoes produce infected offspring, we expect T1h to double. Note that all secondary infections resulting from a single infected vector must occur in the next vector-host-vector transmission cycle; hence, the vector type reproductive number takes the form(34)T1v=R0H+qv.For details on the derivation of both types of reproductive numbers, see Supplementary Materials, Section A-1. A number of methods for deriving R0 exist and in [21] it is demonstrated that the resulting R0 is not unique but they may have a common threshold at R0=1. Hence, their behaviour below and above unity may show significant differences. In Figure 1© we examine how vertical transmission efficiency contributes to the qualitative behaviour of R02,T1h, and T1v above the unity threshold criterion. An important characteristic pattern of their relationship is that below a certain percentage of vertical transmission efficiency the three reproductive numbers are undistinguishable but beyond the relationship between vertical infection and host type reproductive number becomes exponential. This behaviour results from the fact that in each generation, the number of host infections is proportional to the number of infected mosquitoes, which decreases proportionally to the vertical infection efficiency. Hence, for low vertical infection efficiency, possibly the first generations make a significant contribution to subsequent transmission cycles (for more discussion see [23, 27]).In Figure1(b) we compare the basic reproductive number R0 derived through the next-generation method with the one derived from the endemic equilibrium R0e with respect to vertical transmission efficiency. Clearly, the two give different relationships and it is easy to see that the effect of vertical infection is completely diluted on R0. Hence, this measure if used for disease control effort is more likely to provide misleading guidance. The relationship of the basic reproductive number R0e with mosquito to host ratio m and vertical infection qv shows some interesting patterns. For m>1 the relationship between R0e and qv is exponential as described in Figure 1(c); however for m<1 it is linear rendering qv having little effect if not negligible for qv∈(0,0.8) (see Figure 1(d) dashed lines). This suggests that when there are more hosts compared to mosquitoes, vertical infection efficiency becomes completely negligible as in the model settings, and if m<1, sustained host-vector transmissions may be impossible. On the other hand, we observe that for qv⋙0.8 and even for m<1, R0e becomes much larger highlighting the importance of the interplay between vertical infection and the ratio of mosquitoes to hosts.Remark 7.
In the settings of vector-borne diseases both host type reproductive numberT1h or (R0e) and vector type reproductive number T1v are critical for providing useful informative indicator for guiding disease control efforts.


## 3.1.1. Model Equilibria and Stability Analysis

The basic host-vector model with vertical transmission exhibits two equilibria, namely, the disease-freeE10 and the endemic E1∗, respectively. At the disease-free equilibrium,(8)E10=Sh0,Ih0,Iv0=1,0,0,both vector and host populations persist but with no disease. The prevalence of the disease is denoted by(9)E1∗=Sh∗,Ih∗,Iv∗=a+μhR0ea+μhR0h,μhR0e-1cR0e+bμh/V,μhR0e-1a+μhR0e,where a=phvαm, b=pvhα, c=γh+μh, V=μv(1-qv), and(10)R0e=abVc=phvpvhα2mμv1-qvγh+μh.The Jacobian matrix of system (1) at E10 is given by(11)J1′=-μh0-a0-ca0b-V,such that the characteristic polynomial of matrix (11) is then given as(12)P1λ=λ3+a2λ2+a1λ+a0,with a2=V+c+μh, a1=Vc(1-R0e)+(V+c)μh, and a0=μhVc(1-R0e). The coefficient a2>0 and both a1,a0 are nonnegative if and only if R0e<1. Hence, all Routh stability criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Furthermore, for R0e=1, (9) becomes the disease-free equilibrium. Therefore, the following result holds.Theorem 3.
The disease-free equilibriumE10=(1,0,0) exists and it is globally asymptotically stable if R0e⩽1.Alternatively, the global stability of the disease-free equilibriumE10 can be established using the following Lyapunov function:(13)VSh,Ih,Iv=u1Sh-Sh0-Sh0lnShSh0+u2Ih+u3Ivwhere u1,u2,u3 are some positive constants. Calculating the derivative of V along the solutions of system (1), we obtain(14)V′=u1Sh-Sh0ShμhSh0-Sh-aShIv+u2aShIv-γh+μhIh+u3bIh1-Iv-μv1-qvIv,(15)=-u1μhSh0-Sh2Sh-u3μv1-qv1-u1au3μv1-qvIv-u2γh+μh1-u3bu2γh+μhIh(16)-u3bIhIv+u2-u1aShIv.Choosing u1=u2=b/(γh+μh) and u3=1, (16) becomes(17)V′=-u1μhSh0-Sh2Sh-μv1-qv1-R0eIv-bIhIv. Thus, V′ is negative for R0e⩽1. Note also that V′=0 if and only if Sh=Sh0 and Ih=Iv=0. Therefore, the largest invariant set for (1) is the singleton {E10}. Hence, by LaSalle’s invariance principle [34], E10 is globally asymptotically stable when R0e⩽1 and Theorem 3 is valid.Remark 4.
Clearly, the endemic equilibriumE1∗ exists and is unique for R0e>1. This excludes the possibility of occurrence of backward bifurcation. This result is of great epidemiological significance in guiding efforts for disease control as it indicates that R0e=1 is the critical epidemic threshold.To establish the local stability of the endemic equilibrium we evaluate the Jacobian of the system atE1∗, which gives(18)J1∗=-μh+aIv∗0-aSh∗aIv∗-caSh∗0b1-Iv∗-bIh∗-V.The characteristic polynomial of matrix (18) is then given by(19)P2λ=λ3+b2λ2+b1λ+b0,where (20)b2=V2cR0eμh+a+Vc2R0eR0eμh2+a+VaR0eμhvR0e+b+VbR0e2μh+bμh2R0ea+c+μh+abcμh,b1=V2cR0e3μha+c+μh+VbR0e3μh2a+μh+Vc2R0e3μha+μh+VbcR0e2μh2R0e-1+abcR0eμh2R0e-Vμh+bcR0e2μh3,b0=V2c2R0e3μh+VbcR0e2μh2R0e-1+V2ac2R0e2R0e-1+VabcR0eμhR0eR0e-1-1.The coefficient b2>0 and both b1,b0 are nonnegative if and only if R0e>1. Hence, all Routh criteria are satisfied; that is, the three eigenvalues of matrix (11) are negative or have negative real parts. Therefore, the following results holds.Theorem 5.
The endemic equilibriumE1∗=(Sh∗,Ih∗,Iv∗) exists and it is locally asymptotically stable if R0e>1.A global stability result for the endemic equilibriumE1∗ of system (1) is given below.Theorem 6.
IfR0e>1, the endemic equilibrium E1∗=(Sh∗,Ih∗,Iv∗) is globally asymptotically stable.Proof.
LetD1=v1Sh-Sh∗-Sh⁡∗ln⁡(Sh/Sh∗), D2=v2Ih-Ih∗-Ih∗ln⁡(Ih/Ih∗), and D3=v3Iv-Iv∗-Iv∗ln⁡(Iv/Iv∗) be the components of the Lyapunov function(21)USh,Ih,Iv=D1+D2+D3,where v1,v2,v3 are some positive parameters to be chosen later. Differentiating U along the solutions of system (1), we obtain (22)D1′=v1Sh-Sh∗Shμh-μhSh-aShIv=-v1μhSh-Sh∗2Sh+v1Sh-Sh∗ShaSh∗Iv∗-aShIv=-v1μhSh-Sh∗2Sh+Iv∗-ShSh∗Iv-Sh∗ShIv∗-Iv, for v1=1/aSh∗ and μh=μhSh∗-aSh∗Iv∗ at equilibrium. (23)D2′=v2Ih-Ih∗IhaShIv-γh+μhIh=v2Ih-Ih∗IhaShIv-aSh∗Iv∗Ih∗Ih=ShSh∗Iv-IhIh∗Iv∗-Ih∗ShIhSh∗Iv+Iv∗, for v2=v1 and γh+μh=aSh∗Iv∗/Ih∗ at equilibrium. (24)D3′=v3Iv-Iv∗IvbIh1-Iv-μv1-qvIv=bv3Iv-Iv∗IvbIh1-Iv-Ih∗1-Iv∗Iv∗Iv=bv3Iv-Iv∗IvIh-IhIv-IvIh∗1-Iv∗Iv∗,⩽Iv-Iv∗IvIh-IhIv-IvIh∗, for v3=1/b and μv(1-qv)=Ih∗(1-Iv∗)/Iv∗ at equilibrium. Now all together, (25)U′⩽-v1μhSh-Sh∗2Sh-Iv∗IvIv∗+IhIh∗-2-Sh∗ShIv∗-Ih∗ShIhSh∗Iv-Iv∗IvIv∗+IhIh∗-2-2Iv∗+2Iv+Ih-IhIv-IvIh∗-Iv∗IvIh+Iv∗+Iv∗Ih∗,⩽-v1μhSh-Sh∗2Sh-2Iv∗IvIv∗+IhIh∗-2+2Iv1-Iv∗Iv+Ih1-Iv∗Iv+Iv∗Ih∗1-IvIv∗+Iv∗1-Sh∗Sh1+Ih∗Sh2IvIhSh∗2Iv∗. Using the inequality 1-x+ln⁡x⩽0 for x>0 with equality holding if and only if x=1 and the fact that the arithmetic mean is greater than or equal to the geometric mean, we obtain U′(Sh,Ih,Iv)⩽0 for all Sh,Ih,Iv>0. Furthermore, we obtain that U′(Sh,Ih,Iv)=0 holds only when Sh=Sh∗,Ih=Ih∗,Iv=Iv∗ and that E1∗ is the only equilibrium state of these systems on this plane (line). Therefore, by LaSalle’s invariance principle [34], the positive equilibrium E1∗ is globally asymptotically stable.


## 3.1.2. Epidemic Thresholds, Vertical Infection, and Basic Properties

One of the most important critical thresholds in epidemic models is the basic reproductive number,R0, which is usually found using the next-generation method, as the dominant eigenvalue of thenext-generation matrix [35, 36]. Following the method in [36] we write system (1) consisting only of infectious compartments as the difference between new infection and transfer rates and the resulting Jacobian matrices evaluated at the disease-free equilibrium E0=(1,0,0). For system (1) we have two infected classes, namely, Iv and Ih. It follows that the transmission and the transfer matrices F and F, respectively, are defined as(26)F=qvμvpvhαphvαm0,V=μv00γh+μh. Unlike in host-vector models without vertical transmission, the diagonal elements of the transmission matrix F are nonzero. This stems from the fact that in the presence of vertical transmission there is vector to vector transmission, which completely changes the nature of the basic reproductive number. Thus, the next-generation matrix, K, is then given by(27)K=FV-1=qvpvhαγh+μhphvαmμv0, and the resulting dominant eigenvalue of the spectral radius FV-1, which is the basic reproductive number, R0, is given by(28)R0=12qv+12qv2+4R0H,withR0H=pvhphvα2mμvγh+μh.When there is no vertical transmission, qv=0, as in the case of malaria, R0=R0H is simply the geometric mean of the product of the number of new infections in hosts from one infected vector and the number of new infections in vectors from one infected host, in the limiting case that both populations are fully susceptible. The interpretation and epidemiological significance of R0H are well established. It is easy to see that the transmission of infection is increased with efficiency of vector biting α and probabilities of successful infection transmission, but it is hindered by high mosquito death rates and faster host recovery. The biting rate α appears as α2 because it enters twice in the transmission cycle [31]. An important parameter is the ratio of female mosquitoes to hosts m=Nv/Nh, which is central for disease spread according to model settings. R0 increases with the number (or density) of mosquitoes but decreases with the number (or density) of host population. This results from the asymmetry in the dependence of the vector’s biting rate on the sizes of the host and vector populations such that when there are many more hosts compared to mosquitoes, sustained transmission may be impossible. Therefore, in the absence of vertical transmission, for the infection to successfully spread and invade, the ratio of mosquitoes to hosts needs to be sufficiently large so that double bites are common [37]:(29)NvNh>μvγh+μhpvhphvα2, where the critical ratio is given by mc=μv(γh+μh)/pvhphvα2. Note that each mosquito could infect less than one host on average, and yet R0 could still be more than unity. To elucidate this fact we write R0H as a product of each host type single-step reproductive number; that is,(30)R0H=Rhv×Rvh=phvαmμv×pvhαγh+μh, where Rvh represents the number of new infections in mosquitoes from a single infected host while Rhv represents the number of new infections in hosts from a single infected mosquito. Clearly, R0H can be greater than unity even when one of these reproduction numbers is less than unity, and it can also be less than unity even if one of its components is greater than unity. In Figure 1(a) we depict contours plots corresponding to the overall R0 in (31) along the plane (Rhv,Rvh) in the presence of vertical transmission. Its effects in this asymmetric relationship between R0 and its components is not very pronounced and the ratio of mosquitoes to hosts remains one of the leading factors when there is a large disparity between the sizes of the host and vector populations.Figure 1
Relationships among epidemic thresholds with respect to vertical transmission and ratio of mosquitoes to hosts. (a) The relationship between number of new infections in hosts and number of new infections in vectors. (b) Comparison between basic reproductive numbersR0 and R0e derived from two different methods with respect to vertical infection. (c) Comparison of behaviours among the complete transmission cycle R02 and host type and vector type reproductive numbers with respect to vertical transmission efficiency away from the common epidemic criterion. (d) How the ratio of mosquitoes to hosts m and vertical infection qv contribute to patterns of R0e. Dashed lines indicated contours lines for R0e<1 while continuous lines are for R0e⩾1.
(a) (b) (c) (d)Clearly, the geometric mean is less than the average expected number of new infections per generation. This is the case where the next-generation method fails to produce the correctR0 if transmission between hosts is intermediated by another host (for more discussion about the failure of this method see [21]). Instead, it gives the weighted average lying between the number of new infections each individual produces in the next infection event. Note that if the number of hosts is increased, this deficit is compounded. This can have serious implications for guiding disease control efforts as it fails to provide the actual severity of the infection. Therefore, in this study for application purpose we propose the use of another epidemic threshold. This is R0 from (10) derived in Section 3.1.1, which is hereby referred to as the ‘effective’ reproductive number:(31)R0e=11-qvpvhphvα2mμvγh+μh,as (1) it satisfies the property that the endemic equilibrium E1∗ only persists if R0e is greater than unity and (2) the endemic equilibrium exists without occurrence of backward bifurcation, meaning that R0e=1 is the correct critical value. At this point, it is of particular interest to establish the critical ratio of mosquitoes to hosts in the settings of vertical transmission. From (31) we obtain(32)NvNh>1-qvμvγh+μhpvhphvα2, that is, a new critical ratio mc∗=(1-qv)mc<mc for 0<qv<1. This result indicates that, in the presence of vertical transmission, the critical ratio of female mosquitoes to hosts reduces at a rate proportional to vertical transmission efficiency. This indicates an earlier occurrence of the point where an epidemic is just possible as above this level the equilibrium prevalence is expected to rapidly increase to its asymptotic value.Others are epidemic thresholds that provide a direct measure of the control effort required for disease eradication. These measures are known as “type” reproductive numbers and are determined based on the next-generation matrix [38, 39] such that for our case we have host and vector type reproductive numbers T1h and T1v, respectively. Note that in the absence of vertical transmission the host and vector type reproductive numbers are both equal to the square of the basic reproductive number, R0 in (41), which turns out to be R0e in (31) for qv=0. In the presence of vertical infection, R02 gives the expected number of secondary infections after one average, complete (host-vector-host or vector-host-vector) transmission cycle but does not correspond to a specific population type [23]. For the case of host population, one infected host leads to some secondary host infections in the next host-vector-host transmission cycle. This results exactly from the horizontal transmission mechanism and it is given by R0H. However, further secondary host infections may also occur after any number of vector-vector transmission cycles as a result of transovarial transmission. From the expression of R0 in (31) it can be seen that new infections resulting from vector-vector transmission cycle are given by 1/(1-qv); hence the host type reproductive number is then given by(33)T1h=11-qvR0H=R0e.From (33), R0H=R0e if 1/(1-qv)=1; that is, qv∗=0 which is the boundary condition. An important question would be, what is the proportion of vertical infection that doubles R0H? For 1/(1-qv)=2 we obtain qv∗=0.5, which means that if about 50% of infected mosquitoes produce infected offspring, we expect T1h to double. Note that all secondary infections resulting from a single infected vector must occur in the next vector-host-vector transmission cycle; hence, the vector type reproductive number takes the form(34)T1v=R0H+qv.For details on the derivation of both types of reproductive numbers, see Supplementary Materials, Section A-1. A number of methods for deriving R0 exist and in [21] it is demonstrated that the resulting R0 is not unique but they may have a common threshold at R0=1. Hence, their behaviour below and above unity may show significant differences. In Figure 1© we examine how vertical transmission efficiency contributes to the qualitative behaviour of R02,T1h, and T1v above the unity threshold criterion. An important characteristic pattern of their relationship is that below a certain percentage of vertical transmission efficiency the three reproductive numbers are undistinguishable but beyond the relationship between vertical infection and host type reproductive number becomes exponential. This behaviour results from the fact that in each generation, the number of host infections is proportional to the number of infected mosquitoes, which decreases proportionally to the vertical infection efficiency. Hence, for low vertical infection efficiency, possibly the first generations make a significant contribution to subsequent transmission cycles (for more discussion see [23, 27]).In Figure1(b) we compare the basic reproductive number R0 derived through the next-generation method with the one derived from the endemic equilibrium R0e with respect to vertical transmission efficiency. Clearly, the two give different relationships and it is easy to see that the effect of vertical infection is completely diluted on R0. Hence, this measure if used for disease control effort is more likely to provide misleading guidance. The relationship of the basic reproductive number R0e with mosquito to host ratio m and vertical infection qv shows some interesting patterns. For m>1 the relationship between R0e and qv is exponential as described in Figure 1(c); however for m<1 it is linear rendering qv having little effect if not negligible for qv∈(0,0.8) (see Figure 1(d) dashed lines). This suggests that when there are more hosts compared to mosquitoes, vertical infection efficiency becomes completely negligible as in the model settings, and if m<1, sustained host-vector transmissions may be impossible. On the other hand, we observe that for qv⋙0.8 and even for m<1, R0e becomes much larger highlighting the importance of the interplay between vertical infection and the ratio of mosquitoes to hosts.Remark 7.
In the settings of vector-borne diseases both host type reproductive numberT1h or (R0e) and vector type reproductive number T1v are critical for providing useful informative indicator for guiding disease control efforts.


## 3.2. Effects of Ratio of Mosquitoes to Hosts and Vertical Transmission on Disease Prevalence

To fully capture effects of vertical transmission efficiency on vector-borne transmission models, we explore how this parameter influences the disease spread patterns.Results are shown in Figure2. The left figure explores how the first peak of an epidemic denoted by max(Ih(t)) varies with vertical transmission efficiency, qv. The results show that max(Ih(t)) increases linearly with qv. The right-hand figure plots the changes on time of occurrence of each epidemic peak at each value of vertical transmission efficiency. The results suggest that as vertical transmission efficiency increases the time position of each epidemic peak decreases linearly. This suggests that if vertical infection efficiency is increased, the first peak of an epidemic is likely to be attained earlier than the expected time. This suggests that while vertical transmission efficiency increases the epidemic size of an outbreak, it linearly decreases the time of the peak of an outbreak. The lower figures, (c) and (d), describes how the prevalence of infected hosts, Ih, varies with the ratio of female mosquitoes to hosts for different values of vertical transmission efficiency. Figure 2(c) indicates that the prevalence of the disease saturates at larger values of the ratio at earlier stage of the initial epidemic while Figure 2(d) indicates that the prevalence of the disease saturates at lower values of the ratio at later stage of the initial epidemic. This result stems from the fact that at earlier stage of the initial epidemic the size of the vector population is still very high compared to later stages of evolution of the epidemic.Figure 2
The left figure (a) explores effects of vertical transmission on the prevalence of the disease on hosts, where (Ih(qv)) is plotted at max(Ih(t)) for every qv. The right figure (b) explores how the first peak of an epidemic (max(Ih(t))) varies with vertical transmission efficiency (qv). The lower figures represent the dynamics of Ih(t) at time position (78) corresponding to day 27 (figure c) and at time position (150) corresponding to day 50 (figure d) for different values of vertical infection.
(a) (b) (c) (d)


## 3.3. Analysis of the Model with Vector Stages and Vertical Transmission

### 3.3.1. Model Equilibria

Model equilibria are obtained by setting the right-hand side of system (2) equal to zero. The model has three equilibria, two disease-free and one endemic equilibrium. Details about the derivation of the components of the endemic equilibrium are given in Section A-2 of Supplementary Materials. The first trivial equilibrium is both mosquito-free and disease-free:(35)X0=Ls0,Li0,Sv0,Iv0,Sh0,Ih0=0,0,0,0,1,0, which corresponds only to the presence of hosts (this can be livestock for the case of RVF). However, this equilibrium is not biologically feasible or relevant in particular in tropical regions where vector-borne diseases such as Dengue and RVF are endemic. The second equilibrium corresponds to the coexistence of both vectors and hosts but without infection. This is the disease-free equilibrium:(36)X1=Ls1,Li1,Sv1,Iv1,Sh1,Ih1=l∗,0,P∗,0,1,0, which can be used for computation of the basic reproduction number when using the next-generation method [36]. The nonzero components are(37)l∗=δr1-Q0-1δr+θ+μL,P∗=θ1-Q0-1θ+δμv,withQ0=rθμvθ+μL. Hence, X1 is biologically significant if Q0>1; that is, the mosquito population exists and can establish itself if on average each adult female mosquito produces at least one female mosquito at oviposition rate r during her life time period 1/μv after successfully surviving the aquatic stage at probability θ/(θ+μL).The endemic equilibrium which represents the prevalence of the disease in both host and vector populations is given by(38)X∗=Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗, where(39)Ls∗=l∗-Li∗,Sv∗=P∗-Iv∗,Iv∗=μh1-Sh∗λhvSh∗,Sh∗=1-γh+μhμhIh∗,(40)Ih∗=1-qvμvγh+μhR0c-1λvhλhvP∗γh+μh/μh+λvh/1-qvμv,where(41)R0c=λvhλhvP∗1-qvμvγh+μh.Similar to the discussion in Section 3.1.2 the above epidemic threshold in (41) is the “effective” or “control” reproductive number since it satisfies the property that the endemic equilibrium E∗ only persists if R0c is greater than unity and there is no backward bifurcation.Remark 8.
The relationship between the endemic equilibrium andR0 which depends on the vector population threshold indicates that the existence of this equilibrium is governed by the norms of the vector population threshold. Therefore, this vector population threshold is an important parameter in vector control efforts as it provides a means for identifying key factors for reducing the vector population.In the following section we deriveR0 using the next-generation method.


### 3.3.2. The Basic Reproductive Number and Other Control Thresholds

The basic reproductive number,R0, is a concept of great epidemiological significance. Following the method of [36], we write system (2) considering only disease compartments, Li,Iv and Ih as x˙=F-V with xT=(Li,Iv,Ih), where(42)F=δr1-LqvIvθδ1-NvLi+λvhIhSvλhvIvSh,V=θ+μLLiμvIvγh+μhIh. Evaluating the Jacobian of above matrices at the disease-free equilibrium X1=(l∗,0,P∗,0,1,0), the basic reproduction number is then the spectral radius of the next-generation matrix; that is, R0=ρ(F(X1)V(X1)-1),(43)R0=12qv+12qv2+4R0H,withR0H=λvhλhvP∗μvγh+μh.Following the discussion in Section 3.1.2 this epidemic threshold (43) can lead to an underestimate of the control effort required if used to guide disease control strategies. To overcome this inconsistency we use type reproductive numbers to estimate efforts required to provide informative indicators in control of vector-borne diseases. To avoid repetition we omit the steps required in the derivation process, such that the host type reproductive number is(44)T1h=R0H1-qv=R0c,and the vector type reproductive number is then given by(45)T1v=R0H+qv.In the majority of mosquito-borne diseases, treatment or adequate vaccine is not likely to be feasible, either due to financial constraints or nonavailability, in particular in African communities with low resilience to economic challenges. In such cases the main preventive measures are individual protection against mosquito bites and the control of the proliferation of both the larval and adult mosquitoes. This can be achieved by lowering appropriate epidemic thresholds below unity, that is, forcing the reproductive number of the pathogen below 1. In the context of our model this implies reducing either of or both type reproductive numbers (T1h and/or T1v) below unity. Hence, we define the eradication effort to be the percentage reduction in the mosquito population size required to prevent endemic transmission of the pathogen [23], such that if using T1h the eradication effort is 100(1-1/T1h)=100(R0H+qv-1)/R0H and for T1v the effort required to eradicate the disease is 100(1-1/T1v)=100(R0H+qv-1)/(R0H+qv). Therefore, vertical infection at qv% can be responsible for at most qv% of the required eradication effort, provided R0H>1.


### 3.3.3. Stability of the Host-Only System

This fixed pointX0=(Ls0,Li0,Sv0,Iv0,Sh0,Ih0)=(0,0,0,0,1,0) is an equilibrium point of the host-only system given by (46)dShdt=μh1-Sh, which is local and globally asymptotically stable for all t>0 as the derivative μh1-Sh′ is negative.


### 3.3.4. Local Stability of the Disease-Free Equilibrium

The local stability of the disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is established by analysing the eigenvalues of the Jacobian matrix at X1:(47)JE1=-δrP∗-θ+μLδrP∗δr1-l∗δr1-l∗1-qv000-θ+μL0δr1-l∗qv00θδ1-P∗0-θδl∗-μv-θδl∗0-λvhP∗0θδ1-P∗0-μv0λvhP∗000-λhv-μh0000λhv0-γh+μh.Direct computation shows that λ1=-μh is an eigenvalue of matrix (47) and the remaining are solutions of equation(48)λ5+c4λ4+c3λ3+c2λ2+c1λ+c0=0where(49)c4=δrP∗+2θ+μL+γh+μh+θδl∗+2μv>0,c3=1-qvγh+μhμv1-R0c+δrP∗γh+μh+θl∗δ+2μv+θ+μL+μv2+θδl∗2θ+μL+γh+μh+μv+2γh+μhθ+μL+μv+1-qvμvθ+μL-γh+μh>0,c2=1-qvγh+μhμv1-R0cδrP∗+2θ+μL+θl∗δ+μv+1-qvμvθ+μL-γh+μhδrP∗+θ+μL+θl∗δ+μv+γh+μhθ+μL+μv2+δrP∗θ+μLγh+μh+θl∗δ+μv+γh+μhθl∗δ+2μv+μvθl∗δ+μv+θδl∗θ+μLθ+μL+2γh+μh+μv+θδl∗γh+μhμv>0,c1=1-qvγh+μhμv1-R0cδrP∗θ+μL+θl∗δ+μv+θ+μL2+θ+μLμv+2θ+μLθl∗δ+δrP∗θl∗δ+μv+θ+μLθl∗δ1-qvμvθ+μL-γh+μh+γh+μhμv+θ+μLγh+μhθl∗δδrP∗+θ+μL>0,c0=1-qvθ+μLγh+μhμv1-R0cδrP∗θl∗δ+μv+θl∗δμv.By Descartes’ rule of signs, since c4>0, c3>0, c2>0, and c1>0, then all roots of (48) are negative or have a negative real part for c0>0. The coefficient c0 is nonnegative if and only if R0c<1; hence the following result holds.Theorem 9.
The disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is locally asymptotically stable whenever R0c<1.


### 3.3.5. Global Stability of the Disease-Free Equilibrium

To establish the global asymptotic stability (GAS) ofX1 we use results obtained by Kamgang and Sallet [40], which are the extension of some results in [36]. Around the disease-free equilibrium system (2) can be written as(50)x˙S=A1xxS-xE1,S+A12xxIx˙I=A2xxI where xS is the vector representing disease-free compartments (Ls,Sv,Sh) and the vector xI represents the state of infected compartments (Li,Iv,Ih). This requires rewriting system equation around X1, as follows: (51)dShdt=μh-λhvIvSh-μhSh=-μhSh-Sh1-λhvIvSh,dLsdt=δr1-LsSv+1-qvIv-δrLiSv+1-qvIv-θ+μLLs,=-δrSv+θ+μLLs-Ls1+δr1-Ls1Sv-Sv1-δrSvLi+δr1-qv1-LIv.Note that the latter arises from (θ+μL)Ls1=δr(1-Ls1)Sv1. Then, the following matrices are obtained,(52)A1x=-δrSv+θ+μLδr1-Ls10θδ1-Nv-μv000-μh,A12x=-δrSvδr1-qv1-L000-λvhSv0-λvhSv0,(53)A2x=-θ+μLδr1-Lqv0θδ1-Nv-μvλvhSv0λhv-γh+μh. From basic matrix properties and direct algebraic computation it follows that all eigenvalues of matrix A1 are real and negative and both A1 and A2 are Metzler matrices. Hence, system x˙S=A1(x)(xS-xX1,S) is GAS at the disease-free equilibrium xX1,S. To establish the global stability of overall system (50) at X1 conditions of the following theorem in [40] must be satisfied.Theorem 10.
LetΦ⊂U=R+3×R+3. System (50) is of class C1, defined on U if(1)
U is positively invariant relative to (50).(2)
The systemx˙S=A1(x)(xS-xDFE,S) is GAS at xDFE,S.(3)
For anyx∈Φ, matrix A2(x) is Metzler irreducible.(4)
There exists a matrixA¯2, which is an upper bound of the set M=A2(x)∈M3(R)∣x∈Φ¯, with the property that if A¯2∈M, for any x¯∈Φ¯, such that A2(x¯)=A¯2, then x¯∈R3×{0}.(5)
The stability modulus ofA¯2,α(A¯2)=maxλ∈Sp(A2)Re(λ), satisfies α(A¯2)≤0.
Then,X1 is GAS in Φ¯.Proof.
Clearly, conditions(1)-(3) of the theorem have been satisfied. For all x∈Φ, A2(x) is irreducible because (I+A2x)2>0. An upper bound of the set of matrices M, which is the matrix A¯2, is given by matrix A2(x¯), where x¯=(L¯s1,0,S¯v1,0,1,0)∈R3×{0}, with L¯s1=K1 and S¯v1=K2. Similarly matrix A¯2 is irreducible. Recall that from the Perron-Frobenius theorem for an irreducible matrix you get that one of the matrix eigenvalues is positive and greater than or equal to all others, that is, the dominant eigenvalue. Matrix A2 is exactly the matrix used to compute the basic reproductive number, i.e., the dominant eigenvalue. For more details or proof in general settings see [40].Now conditions(1)-(4) have been verified. To check the last condition, we make use of the following Lemma [40].Lemma 11.
LetH be a square Metzler matrix written in block form H=(ABCD), with A and D square matrices. H is Metzler stable if and only if matrices A and D-CA-1B are Metzler stable.MatrixA2(x) in block matrices takes the following components:(54)A=-θ+μL,B=δr1-Lqv0,C=θδ1-Nv0T,D=-μvλvhP∗λhv-γh+μh. Clearly, A is a Metzler stable matrix and(55)D-CA-1B=-μv1-qvλvhP∗λhv-γh+μh is also Metzler stable if μv(1-qv)(γh+μh)-λhvλvhP∗⩾0, that is, 1-R0c⩾0. Therefore, from Theorem 10 and Lemma 11 the following result holds.Theorem 12.
ForQ0>1 there exists a unique disease-free equilibrium X1=(l∗,0,P∗,0,1,0), which is globally asymptotically stable whenever R0c⩽1.Remark 13.
The above result is of great epidemiological importance. It highlights two fundamental indicators:(1) Although vector control is central for control of vector-borne diseases, it does not mean eliminating all the vectors. Note that the effect of this include both infected and noninfected mosquitoes. (2) Disease eradication efforts are independent of the initial sizes of both vector and host populations; however, the ratio between the two populations is a key factor.


### 3.3.6. Local Stability of the Endemic Equilibrium

Results of the local stability of the disease-free equilibriumX1 suggest that for R0c=1 the Jacobian matrix (47) has zero eigenvalue while the remaining eigenvalues are negative or have a negative real part. Since the algebraic computation involved when establishing the stability of the endemic equilibrium through linearization is quite extensive, we employ the centre manifold theory [41]. This theory is used to examine existence of backward or forward bifurcation. The bifurcation occurs at R0c=1, and choosing qv=qv∗ as a bifurcation parameter, then(56)qv∗=1-λhvλvhP∗μvγh+μh,⇒R0H⩽1.The Jacobian matrix J(X1,qv∗) is the same as matrix (47); hence, J(X1,qv∗) has a simple zero eigenvalue λ=0 when R0c=1 and others given as λ=-μh and roots of the polynomial λ4+c4λ3+c3λ2+c2λ+c1=0, where c1,c2,c3,c4 are as defined in Section 3.3.4. Let ω=(ω1,ω2,ω3,ω4,ω5,ω6) be the right eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are derived by solving J(X1,qv∗)×ω=0, which gives(57)ω=F∗ω4,H∗ω4,G∗ω4,ω4,-λhvμhω4,λhvγh+μh where (58)F∗=rθμhQ0-1-qv∗+1-qv∗μhγh+μhθ+μLδrθ+δμvrθ2μh1-Q0-1δr+θ+μL,G∗=1-qv∗μvθθ+μLγh+μhδr+θ+μL-μhθ+μLQ0v-1+μhμvθ+μL2θ+δμv1-Q0qv∗rθθ+δμv1-Q0-1,H∗=δrqv∗θ+δμvθδr+θ+μL,Q0v=rθ1-qv∗μvθ+μL.Clearly, F∗ and G∗ are nonnegative if and only if qv∗⩽1/Q0, Q0v>1, and Q0>1. Note that Q0v gives the average number of emerged infected female mosquitoes produced by one infected female mosquito. Hence, for Q0v>1 infected mosquitoes will persist in the population, leading to disease persistence. Let us denote υ=(υ1,υ2,υ3,υ4,υ5,υ6) as the left eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are obtained when solving [J(X1,qv∗)]T×ω=0, (T represents transpose), which gives(59)υ=0,M∗υ4,0,υ4,0,λvhP∗γh+μhυ4 where M∗=θμv(δr+θ+μL)/δr(θ+μL)(θ+δμv).Computation of Bifurcation Parameters a ∗ and b ∗. Because υ1=υ3=υ5=0 the associated nonzero partial derivatives of f1,f3,f5 are not considered and the remaining nonzero partial derivatives of system f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) are given by ∂2f2/∂x1∂x4=∂2f2/∂x2∂x4=-δrqv∗, ∂2f4/∂x2∂x4=∂2f4/∂x3∂x2=-θ/δ, ∂2f4/∂x3∂x6=λvh, and ∂2f6/∂x4∂x5=λhv.Thus,a∗=2υ2ω1ω4(∂2f2/∂x1∂x4)+2υ2ω2ω4(∂2f2/∂x2∂x4)+2υ4ω2ω4(∂2f4/∂x2∂x4)+2υ4ω3ω2(∂2f4/∂x3∂x2)+2υ4ω3ω6(∂2f4/∂x3∂x6)+2υ6ω4ω5(∂2f6/∂x4∂x5), which after some algebraic simplification gives(60)a∗=-2υ4ω42δrqv∗F∗M∗+δrqv∗H∗M∗+θδH∗+G∗H∗+1-qv∗μvλvhμh-G∗P∗.Given that the vital dynamics of the host are much more slower λvh/μh⋙G∗/P∗ and using the property that υω=1, then a∗<0 for ω4>0. For b∗ the nonzero partial derivative of f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) is ∂2f2/∂x4∂qv∗=δr(1-l∗), such that b∗=υ2ω4(∂2f2/∂x4∂qv∗)=μvυ2ω4. Using the property that υω=1 we obtain that b∗>0. Since a∗<0 for ω4>0 and b∗>0, the model (2) exhibits a forward bifurcation at R0c=1. Therefore, the following result holds.Theorem 14.
The endemic equilibrium pointX∗=(Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗) exists and is locally asymptotically stable for R0c>1.Remark 15.
The conditionQ0v>1 is an indication that vertical transmission is vital for disease long-term persistence. A particular example is the case of Dengue fever and RVF. Upon failure of this condition, it is much more likely that results of Theorem 14 do not hold. Another important feature of Q0v is that vertical infection is proportional to Q0, and hence vector control is still a viable control strategy even for vector-borne disease with transovarial transmission.The above stability results and those in Theorems9, 10, 12, and 14 can be summarized in Table 1 and in the bifurcation diagram in Figure 3.Table 1
Existence and stability of model equilibria.
X 0 X 1 if (0<qv<1) X ∗ if (0<qv<1) I. Q 0 < 1   &   R 0 c < 1 GS DNE DNE II. Q 0 > 1   &   R 0 c < 1 DNE GS DNE III. Q 0 > 1   &   R 0 c > 1 DNE US LS IV. Q 0 < 1   &   R 0 c > 1 US DNE DNE DNE: does not exist, US: unstable, LS: locally stable, GS: globally stable.Figure 3
The bifurcation diagram in (Q0,R0c) plane. The stability of equilibria in each region is given in Table 1.


## 3.3.1. Model Equilibria

Model equilibria are obtained by setting the right-hand side of system (2) equal to zero. The model has three equilibria, two disease-free and one endemic equilibrium. Details about the derivation of the components of the endemic equilibrium are given in Section A-2 of Supplementary Materials. The first trivial equilibrium is both mosquito-free and disease-free:(35)X0=Ls0,Li0,Sv0,Iv0,Sh0,Ih0=0,0,0,0,1,0, which corresponds only to the presence of hosts (this can be livestock for the case of RVF). However, this equilibrium is not biologically feasible or relevant in particular in tropical regions where vector-borne diseases such as Dengue and RVF are endemic. The second equilibrium corresponds to the coexistence of both vectors and hosts but without infection. This is the disease-free equilibrium:(36)X1=Ls1,Li1,Sv1,Iv1,Sh1,Ih1=l∗,0,P∗,0,1,0, which can be used for computation of the basic reproduction number when using the next-generation method [36]. The nonzero components are(37)l∗=δr1-Q0-1δr+θ+μL,P∗=θ1-Q0-1θ+δμv,withQ0=rθμvθ+μL. Hence, X1 is biologically significant if Q0>1; that is, the mosquito population exists and can establish itself if on average each adult female mosquito produces at least one female mosquito at oviposition rate r during her life time period 1/μv after successfully surviving the aquatic stage at probability θ/(θ+μL).The endemic equilibrium which represents the prevalence of the disease in both host and vector populations is given by(38)X∗=Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗, where(39)Ls∗=l∗-Li∗,Sv∗=P∗-Iv∗,Iv∗=μh1-Sh∗λhvSh∗,Sh∗=1-γh+μhμhIh∗,(40)Ih∗=1-qvμvγh+μhR0c-1λvhλhvP∗γh+μh/μh+λvh/1-qvμv,where(41)R0c=λvhλhvP∗1-qvμvγh+μh.Similar to the discussion in Section 3.1.2 the above epidemic threshold in (41) is the “effective” or “control” reproductive number since it satisfies the property that the endemic equilibrium E∗ only persists if R0c is greater than unity and there is no backward bifurcation.Remark 8.
The relationship between the endemic equilibrium andR0 which depends on the vector population threshold indicates that the existence of this equilibrium is governed by the norms of the vector population threshold. Therefore, this vector population threshold is an important parameter in vector control efforts as it provides a means for identifying key factors for reducing the vector population.In the following section we deriveR0 using the next-generation method.


## 3.3.2. The Basic Reproductive Number and Other Control Thresholds

The basic reproductive number,R0, is a concept of great epidemiological significance. Following the method of [36], we write system (2) considering only disease compartments, Li,Iv and Ih as x˙=F-V with xT=(Li,Iv,Ih), where(42)F=δr1-LqvIvθδ1-NvLi+λvhIhSvλhvIvSh,V=θ+μLLiμvIvγh+μhIh. Evaluating the Jacobian of above matrices at the disease-free equilibrium X1=(l∗,0,P∗,0,1,0), the basic reproduction number is then the spectral radius of the next-generation matrix; that is, R0=ρ(F(X1)V(X1)-1),(43)R0=12qv+12qv2+4R0H,withR0H=λvhλhvP∗μvγh+μh.Following the discussion in Section 3.1.2 this epidemic threshold (43) can lead to an underestimate of the control effort required if used to guide disease control strategies. To overcome this inconsistency we use type reproductive numbers to estimate efforts required to provide informative indicators in control of vector-borne diseases. To avoid repetition we omit the steps required in the derivation process, such that the host type reproductive number is(44)T1h=R0H1-qv=R0c,and the vector type reproductive number is then given by(45)T1v=R0H+qv.In the majority of mosquito-borne diseases, treatment or adequate vaccine is not likely to be feasible, either due to financial constraints or nonavailability, in particular in African communities with low resilience to economic challenges. In such cases the main preventive measures are individual protection against mosquito bites and the control of the proliferation of both the larval and adult mosquitoes. This can be achieved by lowering appropriate epidemic thresholds below unity, that is, forcing the reproductive number of the pathogen below 1. In the context of our model this implies reducing either of or both type reproductive numbers (T1h and/or T1v) below unity. Hence, we define the eradication effort to be the percentage reduction in the mosquito population size required to prevent endemic transmission of the pathogen [23], such that if using T1h the eradication effort is 100(1-1/T1h)=100(R0H+qv-1)/R0H and for T1v the effort required to eradicate the disease is 100(1-1/T1v)=100(R0H+qv-1)/(R0H+qv). Therefore, vertical infection at qv% can be responsible for at most qv% of the required eradication effort, provided R0H>1.


## 3.3.3. Stability of the Host-Only System

This fixed pointX0=(Ls0,Li0,Sv0,Iv0,Sh0,Ih0)=(0,0,0,0,1,0) is an equilibrium point of the host-only system given by (46)dShdt=μh1-Sh, which is local and globally asymptotically stable for all t>0 as the derivative μh1-Sh′ is negative.


## 3.3.4. Local Stability of the Disease-Free Equilibrium

The local stability of the disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is established by analysing the eigenvalues of the Jacobian matrix at X1:(47)JE1=-δrP∗-θ+μLδrP∗δr1-l∗δr1-l∗1-qv000-θ+μL0δr1-l∗qv00θδ1-P∗0-θδl∗-μv-θδl∗0-λvhP∗0θδ1-P∗0-μv0λvhP∗000-λhv-μh0000λhv0-γh+μh.Direct computation shows that λ1=-μh is an eigenvalue of matrix (47) and the remaining are solutions of equation(48)λ5+c4λ4+c3λ3+c2λ2+c1λ+c0=0where(49)c4=δrP∗+2θ+μL+γh+μh+θδl∗+2μv>0,c3=1-qvγh+μhμv1-R0c+δrP∗γh+μh+θl∗δ+2μv+θ+μL+μv2+θδl∗2θ+μL+γh+μh+μv+2γh+μhθ+μL+μv+1-qvμvθ+μL-γh+μh>0,c2=1-qvγh+μhμv1-R0cδrP∗+2θ+μL+θl∗δ+μv+1-qvμvθ+μL-γh+μhδrP∗+θ+μL+θl∗δ+μv+γh+μhθ+μL+μv2+δrP∗θ+μLγh+μh+θl∗δ+μv+γh+μhθl∗δ+2μv+μvθl∗δ+μv+θδl∗θ+μLθ+μL+2γh+μh+μv+θδl∗γh+μhμv>0,c1=1-qvγh+μhμv1-R0cδrP∗θ+μL+θl∗δ+μv+θ+μL2+θ+μLμv+2θ+μLθl∗δ+δrP∗θl∗δ+μv+θ+μLθl∗δ1-qvμvθ+μL-γh+μh+γh+μhμv+θ+μLγh+μhθl∗δδrP∗+θ+μL>0,c0=1-qvθ+μLγh+μhμv1-R0cδrP∗θl∗δ+μv+θl∗δμv.By Descartes’ rule of signs, since c4>0, c3>0, c2>0, and c1>0, then all roots of (48) are negative or have a negative real part for c0>0. The coefficient c0 is nonnegative if and only if R0c<1; hence the following result holds.Theorem 9.
The disease-free equilibrium pointX1=(l∗,0,P∗,0,1,0) is locally asymptotically stable whenever R0c<1.


## 3.3.5. Global Stability of the Disease-Free Equilibrium

To establish the global asymptotic stability (GAS) ofX1 we use results obtained by Kamgang and Sallet [40], which are the extension of some results in [36]. Around the disease-free equilibrium system (2) can be written as(50)x˙S=A1xxS-xE1,S+A12xxIx˙I=A2xxI where xS is the vector representing disease-free compartments (Ls,Sv,Sh) and the vector xI represents the state of infected compartments (Li,Iv,Ih). This requires rewriting system equation around X1, as follows: (51)dShdt=μh-λhvIvSh-μhSh=-μhSh-Sh1-λhvIvSh,dLsdt=δr1-LsSv+1-qvIv-δrLiSv+1-qvIv-θ+μLLs,=-δrSv+θ+μLLs-Ls1+δr1-Ls1Sv-Sv1-δrSvLi+δr1-qv1-LIv.Note that the latter arises from (θ+μL)Ls1=δr(1-Ls1)Sv1. Then, the following matrices are obtained,(52)A1x=-δrSv+θ+μLδr1-Ls10θδ1-Nv-μv000-μh,A12x=-δrSvδr1-qv1-L000-λvhSv0-λvhSv0,(53)A2x=-θ+μLδr1-Lqv0θδ1-Nv-μvλvhSv0λhv-γh+μh. From basic matrix properties and direct algebraic computation it follows that all eigenvalues of matrix A1 are real and negative and both A1 and A2 are Metzler matrices. Hence, system x˙S=A1(x)(xS-xX1,S) is GAS at the disease-free equilibrium xX1,S. To establish the global stability of overall system (50) at X1 conditions of the following theorem in [40] must be satisfied.Theorem 10.
LetΦ⊂U=R+3×R+3. System (50) is of class C1, defined on U if(1)
U is positively invariant relative to (50).(2)
The systemx˙S=A1(x)(xS-xDFE,S) is GAS at xDFE,S.(3)
For anyx∈Φ, matrix A2(x) is Metzler irreducible.(4)
There exists a matrixA¯2, which is an upper bound of the set M=A2(x)∈M3(R)∣x∈Φ¯, with the property that if A¯2∈M, for any x¯∈Φ¯, such that A2(x¯)=A¯2, then x¯∈R3×{0}.(5)
The stability modulus ofA¯2,α(A¯2)=maxλ∈Sp(A2)Re(λ), satisfies α(A¯2)≤0.
Then,X1 is GAS in Φ¯.Proof.
Clearly, conditions(1)-(3) of the theorem have been satisfied. For all x∈Φ, A2(x) is irreducible because (I+A2x)2>0. An upper bound of the set of matrices M, which is the matrix A¯2, is given by matrix A2(x¯), where x¯=(L¯s1,0,S¯v1,0,1,0)∈R3×{0}, with L¯s1=K1 and S¯v1=K2. Similarly matrix A¯2 is irreducible. Recall that from the Perron-Frobenius theorem for an irreducible matrix you get that one of the matrix eigenvalues is positive and greater than or equal to all others, that is, the dominant eigenvalue. Matrix A2 is exactly the matrix used to compute the basic reproductive number, i.e., the dominant eigenvalue. For more details or proof in general settings see [40].Now conditions(1)-(4) have been verified. To check the last condition, we make use of the following Lemma [40].Lemma 11.
LetH be a square Metzler matrix written in block form H=(ABCD), with A and D square matrices. H is Metzler stable if and only if matrices A and D-CA-1B are Metzler stable.MatrixA2(x) in block matrices takes the following components:(54)A=-θ+μL,B=δr1-Lqv0,C=θδ1-Nv0T,D=-μvλvhP∗λhv-γh+μh. Clearly, A is a Metzler stable matrix and(55)D-CA-1B=-μv1-qvλvhP∗λhv-γh+μh is also Metzler stable if μv(1-qv)(γh+μh)-λhvλvhP∗⩾0, that is, 1-R0c⩾0. Therefore, from Theorem 10 and Lemma 11 the following result holds.Theorem 12.
ForQ0>1 there exists a unique disease-free equilibrium X1=(l∗,0,P∗,0,1,0), which is globally asymptotically stable whenever R0c⩽1.Remark 13.
The above result is of great epidemiological importance. It highlights two fundamental indicators:(1) Although vector control is central for control of vector-borne diseases, it does not mean eliminating all the vectors. Note that the effect of this include both infected and noninfected mosquitoes. (2) Disease eradication efforts are independent of the initial sizes of both vector and host populations; however, the ratio between the two populations is a key factor.


## 3.3.6. Local Stability of the Endemic Equilibrium

Results of the local stability of the disease-free equilibriumX1 suggest that for R0c=1 the Jacobian matrix (47) has zero eigenvalue while the remaining eigenvalues are negative or have a negative real part. Since the algebraic computation involved when establishing the stability of the endemic equilibrium through linearization is quite extensive, we employ the centre manifold theory [41]. This theory is used to examine existence of backward or forward bifurcation. The bifurcation occurs at R0c=1, and choosing qv=qv∗ as a bifurcation parameter, then(56)qv∗=1-λhvλvhP∗μvγh+μh,⇒R0H⩽1.The Jacobian matrix J(X1,qv∗) is the same as matrix (47); hence, J(X1,qv∗) has a simple zero eigenvalue λ=0 when R0c=1 and others given as λ=-μh and roots of the polynomial λ4+c4λ3+c3λ2+c2λ+c1=0, where c1,c2,c3,c4 are as defined in Section 3.3.4. Let ω=(ω1,ω2,ω3,ω4,ω5,ω6) be the right eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are derived by solving J(X1,qv∗)×ω=0, which gives(57)ω=F∗ω4,H∗ω4,G∗ω4,ω4,-λhvμhω4,λhvγh+μh where (58)F∗=rθμhQ0-1-qv∗+1-qv∗μhγh+μhθ+μLδrθ+δμvrθ2μh1-Q0-1δr+θ+μL,G∗=1-qv∗μvθθ+μLγh+μhδr+θ+μL-μhθ+μLQ0v-1+μhμvθ+μL2θ+δμv1-Q0qv∗rθθ+δμv1-Q0-1,H∗=δrqv∗θ+δμvθδr+θ+μL,Q0v=rθ1-qv∗μvθ+μL.Clearly, F∗ and G∗ are nonnegative if and only if qv∗⩽1/Q0, Q0v>1, and Q0>1. Note that Q0v gives the average number of emerged infected female mosquitoes produced by one infected female mosquito. Hence, for Q0v>1 infected mosquitoes will persist in the population, leading to disease persistence. Let us denote υ=(υ1,υ2,υ3,υ4,υ5,υ6) as the left eigenvector associated with zero eigenvalue of the Jacobian matrix J(X1,qv∗). Its components are obtained when solving [J(X1,qv∗)]T×ω=0, (T represents transpose), which gives(59)υ=0,M∗υ4,0,υ4,0,λvhP∗γh+μhυ4 where M∗=θμv(δr+θ+μL)/δr(θ+μL)(θ+δμv).Computation of Bifurcation Parameters a ∗ and b ∗. Because υ1=υ3=υ5=0 the associated nonzero partial derivatives of f1,f3,f5 are not considered and the remaining nonzero partial derivatives of system f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) are given by ∂2f2/∂x1∂x4=∂2f2/∂x2∂x4=-δrqv∗, ∂2f4/∂x2∂x4=∂2f4/∂x3∂x2=-θ/δ, ∂2f4/∂x3∂x6=λvh, and ∂2f6/∂x4∂x5=λhv.Thus,a∗=2υ2ω1ω4(∂2f2/∂x1∂x4)+2υ2ω2ω4(∂2f2/∂x2∂x4)+2υ4ω2ω4(∂2f4/∂x2∂x4)+2υ4ω3ω2(∂2f4/∂x3∂x2)+2υ4ω3ω6(∂2f4/∂x3∂x6)+2υ6ω4ω5(∂2f6/∂x4∂x5), which after some algebraic simplification gives(60)a∗=-2υ4ω42δrqv∗F∗M∗+δrqv∗H∗M∗+θδH∗+G∗H∗+1-qv∗μvλvhμh-G∗P∗.Given that the vital dynamics of the host are much more slower λvh/μh⋙G∗/P∗ and using the property that υω=1, then a∗<0 for ω4>0. For b∗ the nonzero partial derivative of f=(f1,f2,f3,f4,f5,f6) at (X1,qv∗) is ∂2f2/∂x4∂qv∗=δr(1-l∗), such that b∗=υ2ω4(∂2f2/∂x4∂qv∗)=μvυ2ω4. Using the property that υω=1 we obtain that b∗>0. Since a∗<0 for ω4>0 and b∗>0, the model (2) exhibits a forward bifurcation at R0c=1. Therefore, the following result holds.Theorem 14.
The endemic equilibrium pointX∗=(Ls∗,Li∗,Sv∗,Iv∗,Sh∗,Ih∗) exists and is locally asymptotically stable for R0c>1.Remark 15.
The conditionQ0v>1 is an indication that vertical transmission is vital for disease long-term persistence. A particular example is the case of Dengue fever and RVF. Upon failure of this condition, it is much more likely that results of Theorem 14 do not hold. Another important feature of Q0v is that vertical infection is proportional to Q0, and hence vector control is still a viable control strategy even for vector-borne disease with transovarial transmission.The above stability results and those in Theorems9, 10, 12, and 14 can be summarized in Table 1 and in the bifurcation diagram in Figure 3.Table 1
Existence and stability of model equilibria.
X 0 X 1 if (0<qv<1) X ∗ if (0<qv<1) I. Q 0 < 1   &   R 0 c < 1 GS DNE DNE II. Q 0 > 1   &   R 0 c < 1 DNE GS DNE III. Q 0 > 1   &   R 0 c > 1 DNE US LS IV. Q 0 < 1   &   R 0 c > 1 US DNE DNE DNE: does not exist, US: unstable, LS: locally stable, GS: globally stable.Figure 3
The bifurcation diagram in (Q0,R0c) plane. The stability of equilibria in each region is given in Table 1.


## 3.4. Numerical Results of the Model with Vector Stages and Vertical Transmission

Numerical analysis using reasonable parameter values for RVF is carried out. In Figure4 we plot time series of infected larval, adult mosquitoes and hosts, and their respective phase portraits. The latter are plotted for t∈[050] years for observing long-term dynamics. The parameter values used give R0c>1 and existence of a stable endemic equilibrium. These results not only confirm findings by analytical analysis (Theorems 9 and 12) but suggest that this equilibrium is also global asymptotic stable for the chosen parameter values. This also implies that the endemic equilibrium is stable not only close to the bifurcation point but also for R0c beyond unity (see inner figures, Figures 4(a)–4(c)).Figure 4
Time series and phase diagrams plot. The inner figures describe the temporal evolution ofLi,Iv,Ih after the transient period for R0c=2.1002, while for the outer figures R0c=1.0012. Parameter values used for R0c=1.0012 are μv=1/25,μh=1/(8∗360),phv=0.14,pvh=0.35,α=256/360, γh=1/4,m=1.5,qv=0.2,θ=1/10,μL=1/18,r=1/12,δ=0.4 and parameter values used for R0c=2.1002 are given in Table 1 of Section A-3 of Supplementary Materials.
(a) (b) (c) (d) (e) (f)


## 4. Discussion and Conclusion

One of the main contributions of mathematical epidemic models is that they enhance our understanding of disease transmission and public health planning. This is possible because with models we are able to derive critical epidemic thresholds in terms of model parameters which carry important disease features and key players. One of these measures is the basic reproduction numberR0, which if greater than unity implies that the disease will persist in the population and fade out otherwise. However, a major challenge is to derive the correct R0 for disease models in which infection transmission is intermediated by another host. This is the case of vector-borne diseases, where practical methods such as the next-generation matrix for deriving R0 fails to give the correct value [21]. Rather, it gives the geometric mean of secondary infection per generation. Using the basic model we discussed alternative measures, known as host and vector type reproductive numbers (T1h and T1v, respectively) [38, 39], to be used for guiding disease control efforts. Further, for both the basic and the complex model we showed that the basic reproductive number derived from the existence of endemic equilibrium is the correct epidemic measure and it is equivalent to the host type reproductive number T1h. In fact, this is an alternative method for finding R0 when it is possible to derive the explicit components of the endemic equilibrium point [21, 42]. Then we established important relationships between the reproductive numbers R02,T1h,T1v and vertical infection efficiency. It follows that there is a critical threshold of vertical infection from which the contribution of this transmission route becomes significant. Otherwise, if vertical infection efficiency is very low, its contribution to disease persistence may be insignificant. The work in [23, 27] discusses this feature in more detail and their findings suggest that these results stem from the fact that for low vertical infection efficiency the virus is rapidly lost unless there is a regular disease amplification in the host population. Recent entomological studies have confirmed this situation in particular for RVF where there has been more and more evidence of disease interepidemic activities [7, 17–20], underlined by sporadic outbreaks on livestock at low levels. During this period vertical infection plays a key role for providing protection against chance extinction, but only if it is efficient enough to span several mosquito generations [23], since at the end of unfavourable season the horizontally transmitted epidemic is declining.Also, our results demonstrated that while the relationship betweenR0 and vertical infection efficiency is linear, the relationship between the correct threshold R0e and vertical infection efficiency is exponential. This shows that the failure of R0 as a critical epidemic threshold compounds itself when varying vertical transmission. Therefore, more caution is needed when deriving the basic reproductive number in vector-borne disease systems that include vertical transmission as one of the transmission modes. Another important property to highlight in our models is the female mosquitoes to host ratio which is more likely to be responsible for the asymmetric relationship between the host-to-vector and vector-to-host reproductive numbers as discussed in other studies [28, 43]. We found that in the presence of vertical infection the critical ratio of female mosquitoes to hosts decreases linearly when increasing vertical transmission efficiency. This result indicates an earlier occurrence of the point where an epidemic is just possible as above this level the equilibrium prevalence is expected to rapidly increase to its asymptotic value [37]. Further we analysed how vertical infection affects the prevalence of the disease at the first peak of an epidemic. The results showed that the prevalence at the first peak is positively correlated to vertical transmission efficiency while the time position of occurrence of each peak is negatively correlated to vertical transmission efficiency. These findings suggest that while vertical transmission efficiency increases the epidemic size of an outbreak, it reduces the duration of an outbreak. Additionally, the relationship between the prevalence of the disease and ratio of female mosquitoes to hosts was accessed for different values of vertical infection. Analysis showed that there is an optimal value of the ratio above which the prevalence of the disease decreases. However, such optimal value changes according to the stage of epidemic development, such that at early stage of the epidemic the optimal value is attained much later than that at later stage of the epidemic. Another interesting result is that after the optimal value of the ratio the effect of vertical transmission changes where higher values have less effect compared to lower values of vertical transmission efficiency.Extending the basic model to include larval compartments with maximal capacities for both larval and adult populations we derived important critical thresholds for both vector population and disease system. We found that the model exhibits three boundary equilibria, namely, a mosquito-free system, a mosquito-host system but disease-free, and one where the disease persists in both populations. A detailed stability analysis is presented and numerical simulations are conducted using parameter values relevant to Rift Valley fever. The mathematical results are used to interpret the biological implications of the relationships betweenR0, vertical infection, and the ratio of female mosquitoes to hosts for assessing viable control measures. Our analysis shows that the qualitative behaviours of the system are completely determined by three key quantities: Q0, the vector population threshold; Q0v, the average number of female mosquitoes produced by a single infected mosquito; and R0c, the effective reproductive number (see the bifurcation diagram in Figure 3 and Table 1). These results provide important qualitative understanding of the interaction between vertical infection and ratio of female mosquitoes to hosts on the prevalence of the disease.The mosquito population persists if each female adult mosquito produces at least one larval mosquito after surviving the larval stage. Then, if an infected individual is introduced in the community without any form of protection, the disease will invade the system leading to appearance of an endemic equilibrium if the epidemic criterion is satisfied (R0e>1). The biological implication of this result is that the disease can only invade an area if already colonized by the mosquito population as expected. Hence, measures targeting critical features of the life cycle of the vector are viable disease control strategies. Furthermore, we found that the persistence of the disease is subjected to vertical transmission efficiency governed by the epidemic criterion Q0v>1. This implies that for the disease to persist for each larval mosquito produced, 1/(1-qv) should be infected. This explicit formulation of the threshold conditions in terms of parameters governing the infection transmission process and vector life history is of great epidemiological significance as it allows disease control efforts to be targeted at specific disease and vector stages. We have also shown that the disease-free equilibrium is locally asymptotically stable for R0c<1 and globally asymptotically stable if R0c⩽1. The latter result is of outmost importance, because it shows that if at any time, through appropriate interventions (e.g., destruction of breeding sites, use of insecticides and repellents on the host, or even vaccination), we are able to lower R0c below 1, then the disease will disappear. When the vertical transmission route is possible, this requires also lowering the contribution from this transmission mode. That is, lowering Q0v which implies keeping qv very close to zero. Therefore, the design of control programs should take into account the implications of this mode of transmission in the case of vector-borne diseases that include vertical transmission in the vector population.The goal of this study was to use mathematical models for discussing important properties of vector-borne disease systems and their implications in attempts for control strategies. Though we have used simple but realistic models, we believe that they remain valuable tools useful in the context of research and for providing qualitative understanding of complex processes underlying vector-borne disease transmission systems in particular in the context of transovarial transmission in the vector population. In conclusion, our results show that(1) including larval stages in the vector population when modelling vector-borne disease systems has benefits, as it allows determining conditions for colonization in terms of vector life cycle features; (2) the failure of R0 derived from the next-generation method compounds itself in the presence of vertical transmission; (3) host type reproductive number gives the correct R0 useful for guiding disease control strategies; (4) while vertical transmission efficiency is likely to increase the size of an epidemic, it decreases its duration; and (5) in the presence of vertical infection the critical ratio of female mosquitoes to hosts reduces linearly when increasing vertical transmission efficiency.



---

*Source: 2618985-2018-11-01.xml*