Mathematical modeling, analysis and simulation of the spread of Zika with influence of sexual transmission and preventive measures

AbstractThe Zika arbovirus transmitted by the Aedes aegypti mosquitoes has been shown to be capable of infecting humans via two routes: the bites of infected vectors and through sexual contacts involving infected and non-infected persons. There is no treatment and current prevention or mitigating efforts rely on the use of the Centers for Disease Control and Prevention recommendations including the use of insecticide-treated bed nets (ITN) and indoor residual spraying (IRS). In this work, we investigate via a mathematical model, the role of ITN and IRS as methods for limiting the impact of Zika transmission. We introduce a model that builds on classical SEIR epidemiological single outbreak models. We compute the basic and control reproduction numbers and the final epidemic size in the presence of control measures ITN and IRS. We derive a gross estimate for the rate of sexual transmission, during the initial stages of the outbreak, in terms of prior estimates of the basic reproduction number from related a...


Introduction
Zika is a vector-borne disease transmitted by infected mosquitoes, primarily Aedes Aegypti, an effective agent for the transmission of co-circulating diseases, dengue and chikungunya. Zika's growing incidence and geographic reach compelled the World Health Organization (WHO) to declare Zika virus infections a global health threat (Fauci & Morens, 2016;Kuno, Chang, Tsuchiya, Karabatsos, & Cropp, 1998). The situation seems to be unique for a virus vector-borne disease since the Zika virus can be sexually transmitted within humans. Further, the fact that Zika virus infections have been linked to microcephaly, a neurological disease, in new borns of pregnant Zika infected women as well as to the Guillain-Barre Syndrome, which in extreme cases may cause paraplegia, have made this emergent disease a most important public health challenge in the Americas (Van Den Berg et al., 2014).
Vector control measures are of fundamental importance in the fight against Zika. Measures that show considerable promise include the introduction of genetically modified mosquitoes, which include the engineering of resistant mosquitoes to arbovirus infection. Mathematical models have been used by several researchers to study the transmission dynamics (Brauer & Castillo-Chavez, 2001, 2012Ross, 1911) and control of vector-borne diseases (Blayneh, Cao, & Kwon, 2009;Gao et al., 2016). These models simulate the effect of different control strategies including mosquito control, reduction of contact with mosquitoes, avoidance of sexual contact (for Zika), sound environmental management practices and community education. These types of models are being used to understand not only the dynamics of the spread of vector-borne diseases but also to conduct testbed experiments for the evaluation of the effectiveness of intervention/control measures aimed at ameliorating their impact at the population level, or at higher levels of organization and over multiple temporal and spatial scales. A typical framework for many of the mathematical models include employing classical compartmental models involving an ordinary differential equation (ODE) system using a Susceptible-Exposed-Infected-Recovered (SEIR) compartmental structure for the humans interacting with a Susceptible-Exposed-Infected (SEI) compartmental model for the vector (Brauer & Castillo-Chavez, 2001;Chowell et al., 2007;Yakob & Clements, 2013). These models have been successful for understanding diseases such as dengue and chikungunya that are also spread by Aedes aegypti vector. There is, however, evidence from WHO and CDC that there is also a direct human sexual transmission component to the Zika disease (Kindhauser, Allen, Frank, Santhana, & Dye, 2016), and it is, therefore, important to quantify the role of the Zika sexual transmission through enhanced SEIR-SEI mathematical models, which is one of the main contributions of this work. Given what is known about the disease, it is also important to include multiple infectious stages and preventative measures which are also considered in this paper for the first time for the Zika virus.
Over the last two decades, two prominent approaches for controlling vector populations, recommended by WHO and CDC, involve the use of Insecticide-Treated Mosquito Nets (ITN) and Indoor Residual Spraying (IRS). Using ITN can help reduce contacts between mosquitoes and humans at home. Further, mosquitoes that remain within the boundaries of sprayed homes after their meals can die as a result of IRS. This paper explores the impact of using selected preventive measures such as ITN and IRS in controlling or ameliorating the spread of the Zika virus and makes an indirect effort to assess the potential contributions of sexually-transmitted Zika infections.
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. Section 3 carries out the basic analysis of the model identifying the basic and control reproduction numbers. Estimates of the rate of Zika sexual transmission are derived in order to assess the contribution of this mode of transmission. Section 4 focuses on the study of disease dynamics in the presence of intervention measures. Conclusions are presented in Section 5.

Mathematical model
In this work, a model in the spirit of the Ross' malaria model Ross (1911), which describes the flow of dynamics of individuals and vector as a function of their epidemiological status is introduced. Human and vector populations are differentiated by the use of the subscripts h and m, respectively. Various levels of complexity can be considered including two classes of infections in humans, asymptomatic and symptomatic, which are assumed to be equally infectious and of similar duration (period of infectiousness) and two modes of Zika transmission directly through sexual contact and indirectly via an infected vector (for humans) or via an infected human (for vectors).
The model is organized around the following flow diagram (see Figure 1). Disease transmission dynamics within this model that incorporates interventions is given by the following SEIR/SEI differential equations: In a population of N h = S h + E h + I h,a + I h,s + R h humans and N m = S m + E m + I m adult female mosquitoes, the susceptible humans S h are assumed to be bitten by infectious mosquitoes I m . These susceptible individuals S h move to the exposed class E h after acquiring the disease via sexual transmission or via an infected vector. This transmission is being modelled via the addition of terms directly proportional to the respective infected human classes (asymptomatic I h,a and symptomatic I h,s ) involved in sexual transmission and an infection rate proportional to the infected mosquito population. Specifically, the per-capita rate of Zika transmission through biting from mosquitoes to humans is modelled in terms of a product involving the biting rate b and vector transmission, that is, by b mh = bβ mh N h with β mh denoting the infectiousness of mosquitoes to humans. Analogously, the rate of Zika sexual transmission from the infected to susceptible humans is given by the usual product with b h = a h N h involving a h , the sexual transmission rate of Zika. Note that the exposed category models the incubation period before a human becomes infectious, contrasting the asymptomatic I h,a and symptomatic I h,s categories. Members of the exposed class E h move to become either symptomatic infectious or asymptomatic infectious at a human incubation rate of ν h which is intrinsic human latent period. A fraction q of the latent becomes asymptomatic infectious. Members from each of the asymptomatic I h,a (t) and symptomatic I h,s infectious classes recover with a rates of γ h,a and γ h,s , respectively. The individual vectors in each sub-population are assumed to move from the susceptible class S m to the exposed class E m through biting of an infected human. Vectors of the exposed class E m move to become infectious I m with a vector incubation rate ν m . The vector natural mortality is given by μ m and the vector transmission rate from the infected human to the vector is given by b hm = bβ hm N h where β hm is the infectiousness of humans to mosquitoes.
To incorporate preventive measures into the model, the effect of ITN is introduced in the rates of transmission from the susceptible human class to the exposed human class through a parameter measured as a per cent (1 − ITN). Note that when ITN = 1, 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, the nets have no effect and the disease can spread through both vector and sexual transmission. In addition, the recovered category R h is included to account for partial immunity to the vector after recovering from the infection. Note that we do not include movement of these individuals back to the susceptible class as there is minimal evidence of individuals contracting the disease once they recover. As in the human model, we also incorporate preventive measures ITN and indoor spraying with residual insecticides IRS into the vector model. The effect of ITN is included in the rates through the term (1 − ITN). We also introduce parameters for the removal of mosquitoes denoted by h and j associated respectively with ITN and IRS. To account for a wide range of behaviours, one can let the values of ITN and IRS to range from 0 to 1.

Mathematical analysis
The basic reproduction number denotes the number of secondary infections generated by an infected vector or human when the population being considered is composed of primarily susceptible humans and vectors. R 0 determines whether there is an outbreak or not (Aparicio, Capurro, & Castillo-Chavez, 2002;Brauer & Castillo-Chavez, 2001;Chowell et al., 2007;Smith, O'Nions, Schilling, Unni, & Bender, 1981). And so, if R 0 < 1, the infection dies out without generating an outbreak with an outbreak taking place whereas R 0 > 1, that is, the number of cases at the start of an outbreak grows exponentially (Wallinga & Lipsitch, 2007).
In this section, we will present two key theoretical results. First, we derive the basic reproduction number R 0 for our mathematical model (1)-(8) using the Next Generation Matrix approach. We also derive a gross estimate for the rate of sexual transmission, during the initial stages of the outbreak, in terms of prior estimates of the basic reproduction number from related albeit not sexually transmitted arboviral diseases. This estimate for the initial exponential growth rate provides a measure of relative contribution of sexual contact to the transmission cycle when mosquitoes are also present. While, the effects of the sexual transmission cannot be separated from the vector-borne transmission, this result can help to predict the rate of sexual transmission with some prior knowledge of the basic reproduction number of other arboviral diseases such as chikungunya or dengue that are not sexually transmitted.

Derivation of the basic reproduction number
Let us recall that the proposed mathematical model for human-vector interaction includes sub-populations with different infectious states. Therefore, we will employ a general approach called the Next Generation Matrix approach (Castillo-Chavez, Cooke, Huang, & Levin, 1989;Castillo-Chavez, Velasco-Hernandez, & Fridman, 1994;Diekmann, Heesterbeek, & Metz, 1990) to find the basic reproduction number R 0 which is given by the following theorem. Theorem 3.1: The basic reproduction number R 0 is given by where

Proof:
Given that the infectious states: E h , E m , I h,s , I h,a , I m in Equations (1)-(8), we can create a vector F that represents the new infections flowing only into the exposed compartments given by Along with F, we will also consider V which denote the outflow from the infectious compartments in Equations (1)-(8) which is given by where p = h · ITN + j · IRS . Next, we compute the Jacobian F from F given by Using matrices F and V one can then compute the Next Generation Matrix FV −1 . Letting this matrix can be calculated to be Note that (i, j) entry of the Next Generation Matrix FV −1 is the expected number of secondary infections in compartment i produced by individuals initially in compartment j assuming that the environment seen by the individual remains homogeneous for the duration of its infection. Also, matrix FV −1 is non-negative and therefore, has a nonnegative eigenvalue. The basic reproduction number can then be computed as R 0 = ρ FV −1 which is the spectral radius of the matrix. This non-negative eigenvalue is associated with a non-negative eigenvector which represents the distribution of infected individuals that produces the greatest number R 0 of secondary infections per generation. In order to calculate the eigenvalues of FV −1 , we consider the characteristic equation where λ denotes the eigenvalues of the matrix and I 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 Note that Theorem 3.1 yields a general result for the basic reproduction number R 0 corresponding to the human-vector model given by Equations (1)-(8) that include both sexual transmission and vector transmission. In the absence of one of these, the derived R 0 simplifies which is given in the next two corollaries.

Corollary 3.2:
In the absence of sexual transmission (a h = 0), the basic reproduction number is given byR where R 2 0,a represents the average number of secondary cases produced by an infectious asymptomatic during their infectious period whereas R 2 0,s represents the average number of secondary cases produced by an infectious symptomatic during their infectious period. Corollary 3.3: In the absence of vector transmission (b = 0), the basic reproduction number is given byR A biological interpretation of this basic reproduction number is that an exposed member introduced into a population of S 0 susceptibles becomes infective with probability (1 − q), in which case they cause a h γ h,s infections during an infective period of length 1 γ h,s or becomes asymptomatic with probability q, in which case they cause a h γ h,a infections during an asymptomatic period of length 1 γ h,a .
The following two corollaries from Theorem 3.1 presents upper bounds for the sexual transmission rate a h under various conditions. Corollary 3.4: In the case of complete intervention using insecticide treated nets (ITN = 1), the infection dies out (R 0 < 1) if the sexual transmission rate a h satisfies the following relationship: Remark. Note that for q = 0.5 which implies that there is equal probability for exposed human sub-population to become symptomatic or asymptomatic infectious individuals, this result yields an interesting fact. If the rate of sexual transmission is less than the harmonic mean of the symptomatic and asymptomatic infectious recovery rates, then the infection dies out. Corollary 3.5: In the absence of any intervention strategies (ITN = 0, IRS = 0), the infection dies out (R 0 < 1) if the sexual transmission rate a h and the biting rate of vector b satisfy the following relationship:

Estimating initial rate of sexual transmission
In this section, we will prove a theorem that will help us to estimate the rate of sexual transmission a h during initial growth of the vector-borne disease. For new emerging infections, such as Zika, the available information about the transmissibility of a new infectious disease epidemic is likely to be restricted to daily counts of new cases. It is well known that these counts increase exponentially in the initial phase of an epidemic (Wallinga & Lipsitch, 2007). This knowledge can help us to determine the initial exponential growth rate for the Zika model represented by Equations (1) Proof: Let us consider the linearization of the model Equations (1)-(8) about the diseasefree equilibrium: Letting y = N h − S h and z = N m − S m , one can obtain the following linearization: The characteristic equation corresponding to this system given by det (J − λI) where J is the Jacobian matrix associated with the system (15) and I is a 7 × 7 identity matrix. The equation can be computed using This reduces to the following equation: The solutions of the linearized equations are linear combinations of exponential terms where the exponents are given by the roots of this equation. To find the initial exponential growth rate, one needs to therefore solve the following fifth degree equation given by Recall that our goal, however, here is not to solve this equation but to find an estimate for the rate of sexual transmission a h with some prior information on the basic reproduction numberR 0 . For this we also assume that one can obtain an observed value of the exponential rise ρ from initial data. The equation that we are interested in then becomes Solving for a h then yields Recalling Corollary 3.2, we note that that basic reproduction number in the absence of sexual transmission is given byR 2 0 = R 2 0,a + R 2 0,s . Using Equations (12) and (13) we then haveR which can be rewritten as Substituting (17) into (16) then yields a closed formula for the rate of sexual transmission a h interms ofR 0 given by Remark. Note that using the result in Theorem 3.6, one can also solve for the basic reproduction numberR 0 in the absence of vector transmission in Corollary 3.3.

Numerical experiments
In this section we will implement the solution to the governing differential Equations (1)-(8). The parameters used in the model will be chosen as shown in Table 1. For this work, we have used parameters from various references as indicated in the table.

Dynamics of the human and vector populations
We implement the system (1)-(8) in MATLAB using a fourth order Runge-Kutta method for solving system of ODEs. For our simulations, we choose the initial populations to be   increasing values of a h that are denoted by different colours in the graph a h = 0 (blue), a h = 0.2 (green), a h = 0.4 (red), a h = 0.6 (cyan), a h = 0.8 (magenta) and a h = 1 (yellow).
The dynamics of the human sub-populations are illustrated in Figures 2 and 3 for the cases with and without control measures. The corresponding dynamics for the vector sub-populations are illustrated in Figures 4 and 5, respectively. Figure 2 plots the dynamics of each of the human sub-populations with no preventive measures (ITN = 0, IRS = 0). This includes the Susceptible (Panel B), Exposed (Panel C), Symptomatic (Panel D), Asymptomatic (Panel E) and Recovered (Panel F). Here, the basic reproduction numbers can be computed for the values of a h using (9) to be respectively, R 0 = 2.02, 2.53, 3.13, 3.81, 4.55, 5.33 that is shown in Panel A.       Figure 5, there are no exposed or infected vectors because of complete intervention (ITN = 1, IRS = 1). This implies that the disease can only transmit sexually and will not be vector-borne. As a consequence, there will be no exposed or infected states for the vector and therefore, Panels C and D in Figure 5 are empty plots.  In summary, Figure 2 illustrates that without any intervention strategy (ITN = 0, IRS = 0), the increased rate of sexual transmission leads to a faster transition as each population is also impacted by the vector transmission. Also note the dynamics of infected vector population that is impacted by the dynamics of the corresponding infected individual in the human sub-population. However, Figures 3 and 5 illustrates that with a complete intervention strategy (ITN = 1, IRS = 1), the mosquitoes are completely eliminated (see Panels B, C and D in Figure 5) and hence the dynamics of the human population is only dependent on the rate of sexual transmission.

Influence of preventive measures
Next, we considered an experiment that varied the effects of ITN and IRS, ranging from 0 to 1. We considered a total of 121 different combinations by choosing ITN and IRS to each take values from 0 to 1 in steps of 0.1. For each combination, the basic reproductive number R 0 was computed using (9). Figure 6 displays R 0 as a function of ITN and IRS. In general, we see that as the use of preventive measures increases, the reproductive number gets closer to 1, eventually falling below that critical value. Note that as ITN increases, the value of R 0 decreases eventually below 1. This is more evident in Figure 7 (panel on the left) which suggests that a high value of ITN is required to cause the disease to go from endemic to dying out. By graphing R 0 for a closer set of ITN values in the interval [0.8, 1], we notice that the critical value for ITN can be 0.85 (see Figure 7: panel on the right). This implies that it is possible to bring R 0 < 1 and eradicate the disease by only using ITN, but this requires a moderately large value of ITN.
In Figure 8, we consider the influence of IRS in the absence of effects from ITN. The figure illustrates that as the value of IRS increases, R 0 does decrease, however, at a very slow rate. This goes to show that the influence of IRS is not enough alone to get R 0 < 1 and hence causing the disease to die out.

Vector transmission vs sexual tranmission
Finally, we consider the dependence of R 0 as a combination of the sexual transmission rate a h and the vector transmission b parameters. Figure 9 illustrates the case with no intervention on the left (ITN = IRS = 0) and full intervention on the right (ITN = IRS = 1). From Equations (12) and (13), letting ITN = 0 yields a linear relationship between the basic reproduction rate and the rate of sexual transmission, namely, R 0 = Qa h . This is clearly seen in the graph on the right. For any other value of ITN ∈ (0, 1], this relationship becomes nonlinear as seen in the panel on the left.

Conclusion and future work
In this work, we have considered a new mathematical model for Zika that incorporates both effects of sexual transmission as well as vector transmission. These effects have been studied analytically and numerically both with and without inclusion of any preventive interventions including ITN and IRS. We also derive the basic reproduction number R 0 for the proposed system of ODEs that accounts for both sexual as well as vector transmission effects and is a function of the preventive measures ITN and IRS. Through numerical experiments we have been able to get further insight into thresholds for disease extinction that can contribute to crucial knowledge of disease control, elimination and mitigation of the spread of Zika.
We also derive a new estimate for the rate of sexual transmission a h at the onset of the outbreak of the disease in terms of the basic reproduction number calculated in the absence of the diseaseR 0 . Note that the relative contribution of sexual transmission cannot be determined based on the exponential rise in incidence alone as the effects of the sexual transmission cannot be separated from the vector-borne transmission (Towers et al., 2016). However, prior knowledge of other arboviral diseases such as chikungunya or dengue that are not sexually transmitted, can help provide an useful estimate for the rate of sexual transmission a h for the initial spread of Zika.
As future work, we plan to consider the influence of key parameters in the model in particular the biting rate b and the rate of sexual transmission a h in predicting uncertainty of extinction thresholds using Latin Hypercube Sampling or Partial Rank Correlation Coefficient. We also hope to study the influence of effective education campaigns in limiting the spread of the disease.
Finally, many of the important ideas of mathematical models for epidemiology build on the study of malaria (Ross, 1911). This pioneering work predicts the possibility of controlling malaria by diminishing a population of the mosquitoes below a threshold. While, this prediction was borne out of practice, controlling mosquito population is a very challenging problem, especially, if the mosquitoes adapt to pesticides. In fact, recent reports Yakob & Walker (2016) suggest that the Zika virus is becoming more resistant to common pesticides. We hope that the research findings from this work provide a framework to help improve our overall understanding for the disease persistence in order to control the spread of Zika. undergraduate research experience in Summer 2016. She also thanks her faculty mentor at George Mason University for guidance on this work.

Disclosure statement
No potential conflict of interest was reported by the authors.