Computational and Mathematical Methods to Estimate the Basic Reproduction Number and Final Size for Single-Stage and Multistage Progression Disease Models for Zika with Preventative Measures

We present new mathematical models that include the impact of using selected preventative measures such as insecticide treated nets (ITN) in controlling or ameliorating the spread of the Zika virus. For these models, we derive the basic reproduction number and sharp estimates for the final size relation. We first present a single-stage model which is later extended to a new multistage model for Zika that incorporates more realistic incubation stages for both the humans and vectors. For each of these models, we derive a basic reproduction number and a final size relation estimate. We observe that the basic reproduction number for the multistage model converges to expected values for a standard Zika epidemic model with fixed incubation periods in both hosts and vectors. Finally, we also perform several computational experiments to validate the theoretical results obtained in this work and study the influence of various parameters on the models.


Introduction
Every year over one billion people are infected from vectorborne diseases including malaria, Dengue, chikungunya, schistosomiasis, leishmaniasis, Chagas disease, Zika, and many more. These diseases affect urban as well as rural communities but thrive primarily among communities with poor living conditions. These vector-borne diseases also impose a substantial economic burden on families and governments. Hence understanding the spread of vector-borne diseases has been a major priority for many countries.
Over the years, many factors have contributed to the increase in the spread of vector-borne diseases including the ability of the vectors to adapt to new habitats, ability of the vectors to become drug-resistant, rapid human movement, and changes in policies on control measures. Mathematical models have often been employed to quantify such dynamics of the vector-borne diseases. These models are often described as compartmental model with the populations under study divided into compartments and, with appropriate assumptions, the different subpopulations transfer between these compartments. One of the earliest models was the formulation of a simple SIR model to describe the epidemic [1][2][3], where the entire population being studied was divided into a susceptible class ( ) which consisted of the number of individuals who are susceptible to the disease and are not infected at time , an infected class ( ) which consists of infected individuals who are assumed to be infectious and are able to spread the disease by direct contact with susceptibles, and ( ) which denotes the number of infected individuals who have been recovered and cannot spread the disease again. Most vector-borne diseases include an exposed phase ( ) between being infected and becoming infective. For a vector-borne disease this would mean a SEIR model for humans interacting with an SEI model for vectors as the vectors are not expected to recover in these models in the time span these models are solved. Currently, most models in the literature employ a SEIR-SEI single-stage model where the incubation periods are assumed to be fixed for both humans and vectors. Relaxing this assumption one may also employ a stage-progression model or the so-called linear chain trick that has been used for modeling diseases like HIV and Dengue [4][5][6][7] where the incubation may be modeled as the progression of multiple substages for humans and vectors. While these models have been considered for chikungunya and Dengue ( [8,9] and references therein), they have not been extended to Zika that also includes direct transmission until recently [10,11]. The single-stage model introduced in [10] was extended by [11] to include symptomatic and asymptomatic infectious stages for the human population along with effect of preventative measures.
Two mathematical quantities that are often of interest in these compartmental models is the basic reproduction number and the final size [12]. The basic reproduction number denoted by R 0 is defined as the number of secondary disease cases caused by introducing a single infective individual into a wholly susceptible population of both hosts (humans) and vectors (mosquitoes). Typically if R 0 > 1 an epidemic occurs while if R 0 < 1, there will likely be no outbreak. The value of R 0 helps to quantify the level of control intervention necessary to contain an outbreak. For example, in the case of malaria, a mathematical model was introduced in order to show that malaria can be greatly reduced by reducing the mosquito population density below a certain threshold [13]. There are multiple ways to mathematically estimate the reproduction number for vector-borne diseases [12]. However these estimates vary considerably that may be because of different external factors such as severity of disease, the level of public health surveillance, and local climate condition that can possibly affect the number of vectors and many other such external factors [14].
The final size of the epidemic refers to the number of members of the population who are infected over the course of the epidemic. While there are various approaches to obtain the basic reproduction number corresponding to the model being analyzed, there are no exact solutions for obtaining final size relations for vector-borne diseases. Recently, a final size relation for epidemic models of vectorborne diseases (that also included direct transmission) was obtained for an age of infection model that can be applied to Zika [14]. Specifically, this work derived an upper and lower bound for the final size relation. This model was formulated and analyzed considering infectivity depending on age of infection which allowed arbitrary periods of stay in each compartment and also the inclusion of control measures such as treatment, quarantine, or isolation. While this work provides a new insight to understanding the epidemics of vector-transmitted diseases through a final size estimate, the authors are not aware of any other work that establishes similar estimates with an upper bound and lower bound for a traditional SEIR-SEI vector-borne disease model that includes direct transmission.
For most vector-borne diseases such as Zika, there are currently no vaccines available and resistance to drugs is an increasing threat. Hence the CDC and WHO have recommended vector control as one of the essential ways to prevent disease outbreaks. One such intervention that has shown a lot of promise in vector-borne diseases such as malaria includes using insecticide treated bednets (ITN) which has been proved to be simple, efficient, and costeffective [15][16][17]. Using ITN can help reduce contacts between mosquitoes and humans at home by providing a physical barrier. The insecticide used to treat the bed net also repels mosquitoes ("excito-repellency"or "deterrence") thus increasing the personal protection offered by the net [18,19]. Finally, mosquitoes which are not repelled will most probably be killed as they come in contact with the insecticide as they often rest on the bed net after biting. For Zika there is a need to develop mathematical models that can help provide insight into the relation between increased coverage of ITN and the decrease of disease prevalence through a combination of the personal protection given by the repellency of the insecticide and the community protection given by its insecticidal action.
In this paper, we present the following new contributions for enhancing our understanding of the spread of Zika. First, we build on a single-stage model similar to what is considered in [10] and generalize them by including the impact of using selected preventative measures such as ITN in controlling or ameliorating the spread of the Zika virus. For this model, we derive the basic reproduction number and a sharp estimate for the final size relation. The derivation for the latter is a new alternative to the derivation of the age of infection epidemic model [14]. Specifically, we show that our result matches well with the results presented in [10,14] in the absence of any control measures. Next, we expand the single-stage model to a new multistage model for Zika that incorporates more realistic incubation stages for both the humans and vectors. For this model also we derive a basic reproduction number and a final size relation estimate for the first time. We observe that the basic reproduction number for the multistage model converges to expected values for a standard Zika epidemic model with fixed incubation periods in both hosts and vectors. This is because both the singlestage and the multistage models would be included in an age of infection model. The proof for the final size of the multistage model builds on the derivation for the single-stage model developed in the work and the result applies also to diseases that can be transmitted directly as well as through a vector. The work in this paper incorporates the multistage progression in the intrinsic incubation periods and can be extended to include extrinsic incubation periods as well. Finally, we also perform several computational experiments to validate the theoretical results obtained in this work and study the influence of various parameters in the model.
The outline of the paper is as follows. In Section 2, we present the mathematical framework used to study the transmission dynamics and control of the Zika virus during a single outbreak via a single-stage model. We derive the basic reproduction number and a new upper and lower bound estimate for the final size relation for the single-stage model that incorporates preventative measures such as insecticide treated bed nets. Section 3 carries out the basic analysis for an expanded multistage progression model that incorporates more incubation stages for the humans and the vectors. Finally, we present numerical results in Section 4 for both the Computational and Mathematical Methods in Medicine 3 single-stage and multistage progression models considered in the paper. Discussions and future work are presented in Section 5.

A Single-Stage Zika Model
In this work we develop and analyze an epidemic model for the spread of Zika through both a vector transmission and direct transmission via sexual contact. We will consider a constant total human population size ℎ with ℎ ( ) susceptibles, ℎ ( ) exposed, ℎ ( ) infected, and ℎ ( ) recovered.
Let the rate of Zika transmission through biting from mosquito to human be given in terms of the biting rate which corresponds to the number of bites in unit time. The effective mosquito bites in unit time that a susceptible human receives may be defined as the product of the biting rate and the probability ℎ that a bite transmits the infection which may also be referred to as infectiousness of mosquitoes to humans. Of this a fraction / is with an infective mosquito from ( ). Thus the number of new infective humans in unit time is ℎ ℎ ( / ). Defining the contact rate to be ℎ = ℎ / , the number of new infective humans in unit time is ℎ ℎ . The rate of the spread of Zika through direct sexual transmission from the infected human subpopulation to the susceptible human population is given by ℎ = ℎ / ℎ where ℎ is the sexual transmission rate of Zika. This adds to the number of new infective humans to be ℎ ℎ ℎ .
The vectors are assumed to move from the susceptible class ( ) to the exposed class ( ) through biting of an infected human. For the vectors, we consider a constant birth rate contribution of vectors in unit time and a proportional vector death rate in each of the susceptible ( ), exposed ( ), and infected ( ) vector classes, so that the total vector population size is constant. We will also assume that the vectors do not recover from infection and therefore there is no recovered class for the vectors. The total number of contacts by vectors sufficient to transmit infection therefore is ℎ and the corresponding vector transmission rate from the infected human to the vector is given by ℎ = ℎ / ℎ where ℎ is the infectiousness of humans to mosquitoes. Therefore, the number of new infective vectors in unit time is ℎ ℎ . We assume that the members of the exposed class ℎ move to become infectious at a human incubation rate of ] ℎ which is intrinsic human latent period. Members of the infectious human class recover with a rate of ℎ . Also, we let vectors of the exposed class move to become infectious with a vector incubation rate ] .
To incorporate preventative measures into the model, the effects of ITN are introduced in the rates of transmission from the susceptible human class to the exposed human class through a parameter measured as a percent = 1 − ITN. Note that when ITN = 1 (or = 0), the only movement from susceptible human class to the exposed class is through sexual transmission and not through the vector. On the other hand, if ITN = 0 (or = 1), the nets have no effect and the disease can spread through both vector and sexual transmission. As in the human model, we also incorporate preventative measure ITN into the vector model. We also introduce parameter for the removal of mosquitoes denoted by ℎ associated with ITN. To account for a wide range of behaviors, one can let the values of ITN from 0 to 1.
This leads to the following SEIR/SEI model for Zika transmission: where = ℎ ⋅ ITN. Note that the total human population ℎ appears implicitly in the parameters and there are no new births or deaths in the human population. Moreover, we will assume that the Zika epidemic is started by a visitor from outside the vector population . Hence we will assume (0) = and the total population of the mosquitoes and (0) = (0) = 0.

Derivation of the Basic Reproduction Number.
Recall that the basic reproduction number is defined as the number of secondary disease cases caused by introducing a single infective into a wholly susceptible population of both hosts (humans) and vectors (mosquitoes). Since the proposed mathematical model for human-vector interaction includes subpopulations with different susceptibility to infection, we will employ a general approach called the Next Generation Matrix approach [23][24][25] to find the basic reproduction number R 0 which is given by the following theorem.

Theorem 1. The basic reproduction number R 0 is given by
Proof. Given the infectious stages ℎ , , ℎ , in (1)-(7), we can create a vector F that represents the new infections flowing only into the exposed compartments. The components of the vector F are obtained by considering the terms denoting new infections from the susceptible equations (1) and (5) entering the exposed equations (2) and (6) with ℎ = ℎ and = N .
Along with F, we will also consider V which denote the outflow from the infectious compartments in (1)-(7) which is given by Next, we compute the Jacobian from F given by and the Jacobian from V given by Using matrices and one can then compute the Next Generation Matrix −1 given by ) .
Note that ( , ) entry of the Next Generation Matrix −1 is the expected number of secondary infections in compartment produced by individuals initially in compartment assuming that the environment seen by the individual remains homogeneous for the duration of its infection. Also, matrix −1 is nonnegative and therefore has a nonnegative eigenvalue. The basic reproduction number can then be computed as R 0 = ( −1 ) which is the spectral radius of the matrix. This eigenvalue is associated with a nonnegative eigenvector which represent the distribution of infected individuals that produces the greatest number R 0 of secondary infections per generation. In order to calculate the eigenvalues of −1 , we consider the characteristic equation where denotes the eigenvalues of the matrix and represents the Identity matrix. This can be simplified to yield The characteristic polynomial therefore is the following quadratic equation given by The basic reproduction number R 0 corresponds to the dominant eigenvalue given by the root of the quadratic equation ] . (17) Remark 2. Note that the infected human infects mosquitoes at a rate of ℎ ℎ / over an average time 1/ ℎ which produces ℎ ℎ / ℎ infected mosquitoes. Now a fraction ] /(] + + ) proceeds to become infectious. Next, the infected vectors infect humans at a rate of ℎ / ℎ for an average time of 1/( + ), producing ℎ / ℎ ( + ) infected humans per vector. The result is Also, sexual transmission produces ℎ cases in average time 1/ ℎ which then yields the additional reproductive number The basic reproduction number R 0 for system (1)- (7) can be written in terms of the basic reproduction numbers corresponding to the vector transmission R and direct transmission R as Note that Theorem 1 yields a general result for the basic reproduction number R 0 corresponding to the humanvector model given by equations (1)-(7) that include both sexual transmission and vector transmission. In the absence of one of these, the derived R 0 simplifies to physically meaningful mathematical quantities which are given in the next two corollaries.

Corollary 3.
In the absence of sexual transmission ( ℎ = 0), the basic reproduction number R 0 only corresponds to the vector transmission, that is, Corollary 4. In the absence of vector transmission ( = 0), the basic reproduction number R 0 only corresponds to the direct (sexual) transmission given by Note that in Corollary 3 (in the absence of sexual transmission), the next generation approach employed yields a square root in the reproduction number because it views the transition from humans to vector to humans as two generations.

Final Size for the Single-Stage Model.
In this section we will derive a relation between the basic reproduction number corresponding to the model equations (1)-(7) and the size of the epidemic. Note that the final size of the epidemic, the number of members of the population who are infected over the course of the epidemic, is − ℎ (∞) which is often described in terms of the attack rate (1 − ℎ (∞)/ ℎ ). We will first prove a lemma that will be used to derive the final size relation. (1)-(7), the total number of infected vectors depends on the dynamics of the epidemic and the total number of human infections as follows:
Next, we prove the following theorem that provides an upper bound for the final size for system (1)-(7) in terms of a basic reproduction numberR 0 that is closely related to R 0 that was derived in the previous section. (1)-(7) the final size relation can be bounded above as follows:

Theorem 6. For equations
where the basic reproduction numberR 0 is the sum of the sexual transmission reproduction number R and the vector transmission reproduction number R given byR 0 = R + R .
Remark 7. Note that integrating (1) and adding (1), (2), and (3) from 0 to , we get This leads to the form One can use this implicit relation between ℎ ( ), ℎ ( ), and ℎ ( ) to describe the orbit of solutions. In addition since the right hand side of (30) is finite, the left hand side is also finite and this shows ℎ (∞) > 0.
Next, we prove an estimate that provides a lower bound for the final size relation. The proof relies on the assumption that the vector population has a much faster time scale than the host population and therefore the vector population is at a quasi-steady-stage equilibrium, given by solutions to the equations for , , and in (1)-(7) that are constant functions of , but may depend on ℎ ( ), ℎ ( ), and ℎ ( ).

Theorem 8.
Let the vector population be at a quasi-steadystage equilibrium. The final size relation for (1)- (7) can then be bounded below as follows: where R * 0 is given by Proof. To determine the lower bound for log( ℎ (0)/ ℎ (∞)), we will first obtain a minimum for the Susceptible vector population . Since the vector population is assumed to be at a quasi-steady-stage equilibrium, we leṫ( ) = 0 in (5) which yields This can be rewritten to give since ℎ = ℎ / ℎ and ℎ ≤ ℎ . Therefore, Remark 9. Theorems 6 and 8 provide an upper and lower bound estimate for the final size relation, respectively, to yield Note that for = 1 and = 0, we are able to recover similar estimates that were derived for an age of infection epidemic model that included both vector transmission and direct (sexual) transmission [14].
To summarize, we have introduced the following variations of basic reproduction number in the description for the single-stage model: Remark 10. Note that in the derivation of the final size for the single-stage model we have assumed that (∞) − (0) = 0.
In general, however, one can expect (∞) − (0) = < 0 (as we expect the infections to die out). In this case, one can prove the upper bound estimate (30) as before; however, the lower bound estimate (41) is not a sharp estimate. This is not considered in this paper.
Computational and Mathematical Methods in Medicine 7

A Multistage Progression Zika Model
In this section we extend the single-stage Zika epidemic model to a multistage model by incorporating more realistic incubation period distributions. Specifically, we relax the assumption of fixed incubation period by using a stageprogression model or the so-called linear chain trick [4][5][6][7]. This Zika model incorporates incubation periods as the progression in ℎ incubation substages in humans ( ℎ 1 , ℎ 2 , . . . , ℎ ℎ ) and incubation substages in vectors ( 1 , 2 , . . . , ). This idea of a stage-progression model has been used for Dengue [7] which is caused by the same species of mosquito that causes Zika. The reason for introducing incubation substages is to model (for the first time for Zika) the time between a human is infected and the onset of symptoms due to the infection. These periods are important determinants of the temporal dynamics of the ZIKV transmission and are therefore critical for clinical diagnosis, outbreak investigation, implementation of prevention, programming control measures, and mathematical modeling. Under this formulation, the resulting incubation periods follow a gamma distribution with integer parameters ℎ and , respectively. When the rates of progression between substages are given by ℎ ℎ and for the incubation periods, the resulting gamma distribution has means 1/ ℎ and 1/ for the incubation periods, respectively, and the corresponding variances are given by 1/ ℎ 2 ℎ and 1/ 2 , respectively [7]. In this model we incorporate this stageprogression only in the intrinsic incubation period in humans and vectors. The analysis will be similar if one were to also incorporate stage-progression in the infectious period for the humans also.
We then have the following system of nonlinear differential equations describing the dynamics of Zika through a human-vector interaction aṡ In this section we prove a theorem that provides estimates for the final size for the system (48)-(56). First, we will first prove a lemma that will be used to derive the final size relation.

Theorem 12. For (48)-(52) the final size relation satisfies the following upper bound estimate:
where R 0 , the basic reproduction number corresponding to the system, is the sum of the reproduction numbers corresponding to direct (sexual) transmission R and vector transmissionR which are given by Proof. Adding (48)-(51) yields the following: This implies that ℎ ( )+∑ ℎ =1 ℎ ( )+ ℎ ( ) is a positive decreasing function and therefore the limit exists. The derivative of positive decreasing function tends to zero, and this yields that − ℎ ℎ → 0 and since ℎ > 0, this implies that ℎ → 0.
Integrating (74) we get Noting that ℎ (∞) = ℎ (∞) = 0 and ℎ (0)+∑ ℎ =1 ℎ (0)+ ℎ (0) = ℎ , (75) simplifies to the following: Employing (48): Substituting (76) and simplifying yields Integrating the left hand side and since the rate of sexual transmission ℎ = ℎ ℎ we get Substituting (57) from Lemma 11 into (79) gives us Computational and Mathematical Methods in Medicine 9 where we have used the integral mean value theorem as ℎ ≥ 0 on the interval [0, ∞) with min ≤ * ≤ max ≤ . Using (76) we can finally conclude that where we have used the fact that * ≤ . Define the reproduction numbers corresponding to direct (sexual) transmission R and vector transmissionR , respectively, to be we can now get the following upper bound estimate satisfied by the final size given by where R 0 = R +R . This proves the upper bound.

Corollary 13.
As the number of stages → ∞, the basic reproduction numberR corresponding to only vector transmission converges to the basic reproduction numberR * corresponding to a model that incorporates infected mosquitoes experiencing a fixed incubation period 1/] that is followed by an infectious state from which vectors do not recover. In particular, Proof. Using the definition ofR we havê where = ] /( + ).
Next, we prove an estimate that provides a lower bound for the final size relation. As in the single-stage model, the proof relies on the assumption that the vector population has a much faster time scale than the host population and therefore the vector population is at a quasi-steady-stage equilibrium, given by solutions to (48)-(56) that are constant functions of .

Theorem 14. Let the vector population be at a quasisteady-stage equilibrium. The final size relation for equations (48)-(56) can then be bounded below as follows:
where R * 0 is given by Proof. To prove the lower bound, recall that Also, since the vector population is assumed to be at a quasisteady-stage equilibrium, settinġ= 0 in (53) yields since ℎ = ℎ / ℎ and ℎ ≤ ℎ . Substituting (89) in (88) then gives (91) Remark 16. Note that in the derivation of the final size for the multistage model we have assumed that (∞) − (0) = 0.
In general, however, one can expect (∞) − (0) = < 0 (as we expect the infections to die out). In this case, one can prove the upper bound estimate (72) as before; however, the lower bound estimate (86) is not a sharp estimate. This is not considered in this paper.
To summarize, we have introduced the following variations of basic reproduction number in the description for the multistage model:

Computational Experiments
In this section, we validate our theoretical results developed in this work and perform simulations to predict the dynamics and estimate the basic reproduction numbers and final size relations. We implement the solution to the single-stage system (1)-(7) and the multistage system (48)-(56) in MAT-LAB using a fourth order Runge-Kutta method for solving ordinary differential equations. Specifically, the single-stage system of differential equations is solved using the script ode45 from MATLAB [26].
For our simulations, we considered the total human and vector populations to be ℎ = 1000 and = 4000, respectively, as in [14]. We considered one infective human initially and that the Zika epidemic is started by a visiting vector from outside the vector population . For the parameters in the models, we refer mainly to the data used in [10,14] that corresponded to the 2015 Zika outbreak in Barranquilla, Colombia. With the values for the parameters in Table 1, a linear relation was proposed in [14] given by We note that as the value of ℎ increases, the nature of the graphs also change. Specifically note that the final value of the number of susceptible humans ℎ ( ) decreases as ℎ increases.
Next, we consider the effect of the insecticide treated nets with no sexual transmission. So we let ℎ = 0 and let ITN = 0.1, 0.2, 0.3. The results are illustrated in Figures 4, 5, and 6, respectively. We clearly note that as the value of ITN increases in percentage (10%, 20%, 30%) of protection because of the nets, it takes more time for humans to get infected. Clearly this shows the effect of using nets.
Next, we consider the influence on the rate of sexual transmission ℎ on the various basic reproduction numbers that were obtained for the single-stage model. These include R 0 derived using the Next Generation Matrix approach in (20), the basic reproduction number corresponding to purely direct (sexual) transmission R defined in (19), the basic reproduction number corresponding to purely vector transmission R defined in (18), and the basic reproduction    numberR 0 obtained in Theorem 6 which is the sum of R and R . Figure 7 compares these reproduction numbers as ℎ increases from 0 to 0.413. Next, we perform numerical experiments to validate the upper bound and lower bound estimate for the final size relation obtained in (46). This is illustrated in Figure 8 which clearly demonstrates the validity of the estimate. One may also note that log( ℎ (0)/ ℎ (∞)) is closer to the upper bound as pointed out in [14]. Also, for a purely direct (sexual) transmission,R 0 =R * 0 and therefore we get equality in the estimate that is denoted by the convergence at ℎ = 0.413 in Figure 8. Next, we considered the sharpness of the upper bound estimate by solving for ℎ (∞) in Specifically, we used Newton's method for solving nonlinear equations. The comparison of the results obtained from the simulations of the dynamics of (1)- (7) with the results obtained via Newton's method is illustrated in Figure 9. Note that the -axis values are normalized.
Since it is well known that increased coverage of ITN decreases vector prevalence, we performed our next computational experiment to explore the influence of using increasing values of ITN on the dynamics. The range of values included the absence of nets (ITN = 0) to completely protective nets (ITN = 1). The results are illustrated in Figure 10 which plots the final number of susceptible humans ℎ (∞), the basic reproduction numberR 0 , and the attack rate 1 − ℎ (0)/ ℎ (∞) for increasing values of ITN. The value of ℎ was chosen to be 0.2 for this simulation that corresponds to the inclusion of both vector and direct transmission in the model. As expected, Figure 10 illustrates the usefulness of insecticide treated bed nets to control the epidemic. The figure also can potentially help government officials to decide the level of control measures through insecticide treated bed nets that is needed in a certain area to contain the spread of the epidemic.
Next, the effect of employing a multistage model on the value of the basic reproduction number R 0 derived in Theorem 12 is illustrated in Figure 11 for increasing values of ℎ = 0, 0.1, 0.2, 0.3, 0.4, 0.413. As the number of stageprogression increases the value of R 0 converges for different  values of ℎ . Finally, a convergence study was performed for the full multistage system (48)-(56) for increasing number of stages doubling each time, to demonstrate that the multistage model parallels the single-stage model. Figure 12 illustrates this convergence for fixed value of ℎ = 0.2. Note that the -axis values of the human population are normalized.

Discussion and Future Work
Over the last several decades, there has been an explosion in the development of mathematical models for outbreaks of infectious diseases that has become part of assessing epidemiological phenomena and making health policy decisions [1,12,27,28]. The outbreak of any new disease has always provided both an opportunity and a challenge for mathematicians and scientists. The opportunity leads to the development of improved models and the challenge is to make sure that these models represent reality. Zika is a great example of such as disease [29]. While there is a lot of published work and information available on models, methods, and simulations on infectious diseases that are purely vector-transmitted such as malaria, Dengue, and chikungunya and diseases that spread through direct transmission only such as influenza and AIDS, new diseases such as Zika that includes both vector transmission and direct transmission has provided a challenge for new mathematical models. For example, one of the essential mathematical tools in understanding disease dynamics is the calculation of the basic reproduction number that can help make informed decisions on whether there will be an epidemic or not. Another important quantity is the final size relation that provides a useful relation between the basic reproduction number and the size of the epidemic. While there has been a lot of progress made in computing these for various types of diseases, there is still a lack of complete understanding of these quantities for vector-borne diseases such as Zika that also includes direct (sexual) transmission.
Based on the approach we take, one can obtain different measures of the basic reproduction number for vector-borne diseases. However, there are no exact analytical solutions for the final size relation for such diseases. Nevertheless, one can obtain sharp estimates with an upper and lower bound for the epidemic size [14]. The latter formulated and analyzed a model with infectivity depending on age of infection. The work provided a useful upper bound and lower bound estimate for this age of infection model that applies to vector-borne diseases such as Zika that also includes direct transmission. In this work, we provide an alternate approach to determining similar estimates for the final size relation for enhanced SEIR-SEI single-stage and multistage progression models. The multistage model for Zika considered herein includes multiple incubation substates and was motivated by similar models for Dengue that did not account for the direct (sexual) transmission. Towards this end, we are able to successfully derive a new upper bound and lower bound estimate for the final size relation for the models considered. Moreover, we are able to show that the basic reproduction number for the multistage model proposed converges to the basic reproduction number corresponding to an equivalent nonlinear system of delay differential equations with fixed incubation periods in the humans and vectors.
Another contribution in this work is the inclusion of insecticide treated nets (ITN) that offer a mix of personal protection-blocking the bites of mosquitoes, thereby reducing the transmission from mosquitoes to humans-and community protection-reducing the longevity of mosquitoes and therefore the prevalence of the infectious stage of the disease, in mosquitoes. All the results developed in this work including the basic reproduction number and the final size relation estimate incorporate the influence of the ITN which are recommended by both the CDC and WHO as effective control measures.
We hope that the models, methods, and result from this work can help provide more insight into the propagation of a disease like Zika. The work also provides some opportunities for new avenues for future research. The multistage progression model in this work only considers multiple substates for the incubation for humans and vectors. However, one may also extend this work to also incorporate multiple infectious states in humans and compare such models against standard epidemic models with fixed incubation periods in both hosts and vectors and an exponentially distributed infectious period in hosts. One may also consider including indoor residual spraying (IRS) that can provide a coating of the walls and other parts of a house with a residual insecticide that can kill the vectors when they come in contact. One of the assumptions in this work involves employing a removal rate of vectors, and we need to assume vectors are in a quasi-steady state equilibrium depending on in order to continue to assume a balance relation for the host and vector contact rates. The parameter indicates the dependence of the model on bednets; it could also include a rate of killing of vectors by spraying, or we might wish to include separate control parameters for the effect of bed nets and the effect of spraying. Then an interesting question would be how the basic reproduction number and the final size of the epidemic depend on these control parameters. Finally, to model the effect of bednets, one may need some terms in the host equations that could resemble those in vaccination models. For example, this can include two susceptible compartments, with a rate of transfer from susceptibles without bednets to susceptibles with bednets who would have a smaller contact rate. In this case the rate of using bednets would be a control parameter. Finally, the derivations presented herein assumed that the initial population and the steady state population of the mosquito are the same. One can relax this assumption and as pointed earlier one can rederive the estimates. While the upper bound estimate can be derived in a similar fashion as in this paper, obtaining a sharp lower bound estimate will require some work. All these features and extensions will be considered in a forthcoming paper.

Conflicts of Interest
The authors declare that they have no conflicts of interest.