Using mathematical modelling to investigate the effect of the sexual behaviour of asymptomatic individuals and vector control measures on Zika

Zika is a vector borne disease for which the latest world wide outbreak inspired a renewed interest in epidemiological modelling of vector borne diseases. However, due to the possibility of sexual transmission and the high proportion of asymptomatic individuals, models for similar diseases, such as dengue or chikungunya, are no longer applicable. It is of interest to study how the existence and behaviour of asymptomatic individuals and the potential of them transmitting the disease affect the overall epidemic dynamics. The model presented here aims to be as simple as possible, while at the same time taking into account the features that make Zika unique among other vector borne diseases. This model allows for the exploration of sexual transmission and how the sexual behaviour of asymptomatic individuals may affect the spread of the disease. In addition, the model was used to determine the basic reproductive number, with and without the effect of sexual transmission as well as to implement a simple version of control using Wolbachia bacterium.


Introduction
Zika is a vector-borne disease which is primarily transmitted through the bite of an infected female Aedes mosquito, mainly Aedes aegypti (Gao et al., 2016). Aedes aegypti mosquito is the same mosquito that can transmit dengue, chikungunya and yellow fever. Unlike dengue, chikungunya and yellow fever, however, Zika can be transmitted through sexual contact (Allard, Althouse, Hébert-Dufresne, & Scarpino, 2017). The most common symptoms of the Zika virus infection are fever, rash, headache, joint pain, conjunctivitis and muscle pain. Symptoms usually last anywhere from 3 to 14 days (Krow-Lucal, Biggerstaff, & Staples, 2017). Yet not everyone who is infected with Zika displays or experiences any symptoms and, in fact, it is estimated that 80% of people infected with Zika are asymptomatic (Duffy et al., 2009).
Many of the early models of Zika epidemics assumed that sexual transmission of the virus can be neglected, and that we can rely on previously analysed models of other vector borne diseases. Papers such as Maxian, Neufeld, Talis, Childs, and Blackwood (2017) concluded that sexual transmission may play a crucial role in the existence of Zika outbreaks, and offers conditions for which sexual transmission is relevant.
Later models tend to consider sexual transmission but group symptomatic and asymptomatic humans in one class. A few recent papers include both the asymptomatic and symptomatic classes, as well as sexual transmission (Maxian et al., 2017;Padmanabhan, Seshaiyer, & Castillo-Chavez, 2017). As it is estimated that the vast majority of people infected with Zika do not present any symptoms, splitting the infected population into these two classes offers the opportunity to better understand the dynamics of the disease and the effect that the behaviour of asymptomatic individuals will have on the spread of Zika. If an infected individual does not present symptoms, it is more likely that this individual will not take the precautions necessary to avoid spreading the disease. For this reason, we present and analyse a model which includes both sexual transmission, and symptomatic and asymptomatic classes and analyse how the sexual behaviour of asymptomatic individuals affects the spread of Zika.
This is similar to the model presented in Padmanabhan et al. (2017), but without World Health Organization (WHO) and Center for Disease Control (CDC) preventative measures. Instead, the authors here are interested in investigating the effects of a different vector control measure, Wolbachia. Given that current strategies to battle the spread of Zika such as insecticides or larval biological control have been proven unsustainable and ineffective, new approaches are needed to halt the spread of the disease. The Wolbachia bacterium has been shown to effectively prevent the transmission of Zika and has the ability to reduce the mosquito population (Walker et al., 2011). Wolbachia is naturally found in at least 40% of all known terrestrial insect species, and while it is not normally present in A. aegypti, it can be introduced to them (McMeniman et al., 2009). It is also known that Wolbachia does not affect humans (Dutra et al., 2016).
There are two ways Wolbachia can be beneficial in reducing the spread of Zika. The first is that it can be used to reduce the mosquito population as well as mosquito life span (Schraiber et al., 2012;Walker et al., 2011). When a male mosquito is infected with a strain of Wolbachia its sperm becomes modified such that embryos die during early embryonic development. This method is very much dependent on the prevalence of Wolbachia in the mosquito population, since uninfected female mosquitoes are less likely to mate with Wolbachia-infected males. Thus there is a risk of the strain becoming obsolete (Jiggins, 2017). The second way Wolbachia can help in stopping the spread of Zika is that it has been shown to greatly reduce the capacity of Zika-infected mosquitoes to harbour and transmit the disease (Hughes & Britton, 2013;Ndii, Hickson, & Mercer, 2012).
Numerous projects and studies have already been conducted introducing Wolbachiainfected mosquitoes into populations (Debug Project, 2017;Dutra et al., 2016;Hughes & Britton, 2013;Jiggins, 2017;McMeniman et al., 2009;Ndii et al., 2012;Schraiber et al., 2012 among many others). Most studies which include mathematical modelling of epidemic dynamics and which focus on the use of Wolbachia to control the spread of diseases utilize highly complex models highly focused on the mosquito population (Ferguson et al., 2015;Hughes & Britton, 2013;Ndii et al., 2012), but not many models focus on Zika spread among humans and the effect that Wolbachia can have in the spread of the disease.
Mathematical modelling has played a crucial role in understanding the dynamics of an epidemic outbreak and in developing and testing different control measures (Brauer & Castillo-Chavez, 2012). Since the most recent Zika outbreak in 2015, there was an increase in number of mathematical models of Zika (see Agusto, Bewick, & Fagan, 2017;Bates, Hutson, & Rebaza, 2017;Gao et al., 2016;Maxian et al., 2017;Olawoyin & Kribs, 2018;Padmanabhan et al., 2017;Towers et al., 2016;Wang, Zhao, Oliva, & Zhu, 2017 among many others). In addition, a wide variety of different mathematical models of other vector-borne diseases are common in the literature (Brauer & Castillo-Chavez, 2012). We present a simple model with the goal of showing this model is sufficient to assess the effect of this kind of vector control on the spread of the disease. Local sensitivity analysis for the model parameters is done and, in addition we use Latin Hypercube Sampling to estimate parameters for which only ranges are known.
The paper is organized as follows: in Section 2, a mathematical model is presented in which the human population is divided into five classes while the vector population is divided into three classes with all classes and parameters defined in Section 2. In Section 2.1, we take an epidemiological interpretation of the system to obtain the nextgeneration matrix, and thus, calculate the basic reproductive number. Section 2.2 presents the results for the local sensitivity indices of R 0 . Section 3 deals with vector control and the effect of the use of Wolbachia in controlling the spread of Zika. Lastly, Section 5 summarizes what was presented in the paper and explores new and future directions of the project.

Mathematical model
A deterministic model of transmission dynamics is developed for the Zika virus. The model includes vector to human, human to vector contact and human to human unprotected sexual contact. Vertical transmission is neglected in the model as it is almost negligible (Center for Disease Control and Prevention, 2018; Kucharski et al., 2016;Olawoyin & Kribs, 2018). The H and V subscripts in the variables denote human and vector populations, respectively. The human population follows an SEIR compartmental model and the vector population follows an SEI compartmental model. See Table 1 for a list of the variables in the model.
A member of the susceptible class, S H , moves to the exposed class, E H , after being bitten by an infectious vector or after unprotected sexual contact with an infectious symptomatic or asymptomatic human, which are considered separate classes; I HS and I HA , respectively. The rate of transfer from the susceptible class to the exposed class by means of vector to human contact is given by a H b, where a H is the transmissibility probability from mosquito to human and b is the bite rate. The rates at which susceptible humans are exposed to the virus due to unprotected sexual encounters with symptomatic and asymptomatic infectious humans are ψ S and ψ A , respectively. The parameters ψ S and ψ A are the products of the rate of sexual intercourse and the probability of being infected with the virus through an unprotected sexual encounter. It is assumed that symptomatic and asymptomatic humans are equally infectious but that their sexual habits differ as subjects with symptoms might be less likely to be sexually active and more likely to use protection if they are aware that they are infectious. However, because the discovery of sexual transmission is fairly recent, the transmission rate is still unknown. Once exposed to the virus, humans move to the infectious classes after a latency period of β −1 H days. The exposed individuals will then become symptomatic with probability q, and asymptomatic with probability 1−q. It is assumed that both symptomatic and asymptomatic individuals recover at the same rate γ (Agusto et al., 2017;Caminade et al., 2017). Once a human moves to the recovered class, R H , they can no longer infect a mosquito or another human. This life-long immunity, although still not proven, is an assumption made by many other models (Agusto et al., 2017;Caminade et al., 2017;Manore & Hyman, 2016;Momoh & Fügenschuh, 2018;Padmanabhan et al., 2017;Towers et al., 2016). The death rate due to the Zika virus is assumed to be negligible and is therefore not included in the model (Chikaki & Ishikawa, 2009;Momoh & Fügenschuh, 2018). Given that the model simulates a single possible outbreak, the birth and death rates of humans are also neglected (Ding, Tao, & Zhu, 2016;Kucharski et al., 2016;Towers et al., 2016;Yakob & Clements, 2013). Thus the total human population is assumed constant, N H , as is equal to the sum of the populations of humans in all other classes, i.e. N H = S H + E H + I HA + I HS + R H .
The total population of vectors (mosquitoes), N V , is divided into three classes: susceptible, exposed and infectious. Death rates across each class are assumed to be the same. Mosquitoes are assumed to be born only into the susceptible class at a rate of λ V . Susceptible mosquitoes become exposed after biting an infected (symptomatic or asymptomatic) human at a rate of a V b. After a latency period of β −1 V days, mosquitoes in the exposed class transition to the infected class. Given the short lifespan of mosquitoes, we assume infectious vectors remain infectious for the rest of their life.
See Figure 1 for the model diagram and Table 2 for the list of parameters and their units. The complete model is presented in the system of equations (1). Figure 1. Black arrows represent the flow within each population. Dashed arrows represent the contact between humans and vectors. A human in the Susceptible class, S H , will transition to the Exposed class, E H , when they are bitten by an infectious mosquito or when they engage in unprotected sex with infectious individuals. A mosquito in the Susceptible class, S V , can transition to the Exposed class, E V , when they bite an infectious individual. Exposed humans and vectors are not able to transmit the infection but will eventually transition to the infectious class (I H and I V respectively). The birth and death of vectors are incorporated in this timespan. Birth and deaths of human populations are ignored since the model aims to predict the behaviour of the disease during a single season, and therefore, changes in the human population are negligible.

Calculation and analysis of R 0
The basic reproduction number, R 0 , is defined as the expected number of secondary cases produced by introducing a single infected individual into a completely susceptible population during the entire period of infectiousness of such individual. If R 0 > 1, the disease can invade and potentially turn into an outbreak or epidemic. If R 0 < 1, the disease will be unable to invade and will not produce an outbreak or epidemic (Cushing & Yicang, 1994;Diekmann, Heesterbeek, & Metz, 1990).  Trpis and Häusermann (1995) Note: Parameters with * were estimated (by the median) using Latin Hypercube Sampling given the ranges provided by the literature.
The next-generation matrix (NGM), K, introduced in Diekmann et al. (1990) is the basis for the definition and calculation of R 0 : R 0 = ρ(K). The ijth element of K is the expected number of new cases with state-at-infection i, generated by one individual who has just entered state-at-infection j in a fully susceptible population (Diekmann, Heesterbeek, & Roberts, 2009). States-at-infection are those in which an individual can be in immediately after becoming infected. For our model, states-of-infection are E H and E V .
For the linear algebra approach for computing the NGM, see Appendix A.1. However, here we derive the NGM K in this systematic manner by epidemiological reasoning. We refer to state-of-infection E H with index 1 and E V with index 2. Thus, for the element K 11 , we begin with one individual who has just entered state E H and follow the individual for the remainder of their infectious life to determine how many cases of exposed humans they are expected to produce. We are assuming that with probability 1 the individual survives the E H state and moves to one of the two infectious states, I HS or I HA . Thus the individual moves to state I HS with probability q. While in state I HS , the individual is expected to produce new cases at a rate ψ S (N H /N H ), for an expected time 1/γ . Similarly, the individual moves to state I HA with probability 1−q, expected to produce new cases at a rate ψ A , for an expected time 1/γ . Thus For the element K 12 , we begin with a vector that has just entered state E V and follow it for the remainder of their infectious life to determine how many cases of exposed humans they are expected to produce. Before the vector can infect, it has to survive the exposed state and move to the infectious state. This occurs with probability β V /(β V + μ V ). While in the infectious state this vector is expected to produce new cases of exposed humans at a rate a H b, for an expected time 1/μ V . Therefore Similarly to determine K 21 , we see that an exposed human moves to state I HS with probability q, is expected to produce new cases of exposed vectors with rate a V b, for an expected time 1/γ . Including the movement from E H through I HA we obtain Finally, K 22 = 0 since there is no transmission between vectors. Therefore, the NGM is Thus the reproduction number of the system is given by In addition, if sexual transmission is ignored, then tr(K) = 0 and the reproductive number due to vector-borne transmission only, R 0V , is given by If vector transmission is ignored, then det(K) = 0 and the reproductive number due only to sexual transmission, R 0S , is given by the expression In order to provide a more robust estimation of the reproductive number, Latin Hypercube Sampling (LHS) maximin criteria was performed (Blower & Dowlatabadi, 1994;van  den Driessche, 2017). Based on the ranges for the parameters with a * in Table 2, one million randomly generated set of parameters were used to calculate R 0 . The parameters were generated assuming a uniform distribution. The median and the percentiles at 2.5% and 97.5% are provided (see  Table 2).
Given the values of R 0 , R 0V and R 0S , we calculate the ratios 14 which help us conclude that although the largest contribution to R 0 clearly comes from the vector-human transmission, sexual transmission cannot be ignored.

Local sensitivity analysis of R 0 on different parameters
Local sensitivity analysis demonstrates how the basic reproductive number will change in response to changes in the model parameters. The sign of the index indicates the direction of the response. That is, if the sensitivity index for a given parameter is positive, then an increase in the parameter will increase the value of R 0 . The magnitude of the sensitivity index indicates the relative importance each parameter has on the model's predictions (Moreno et al., 2017). The normalized sensitivity index (a.k.a. elasticity index) of a variable (R 0 ) to a parameter is the ratio of the relative change in the variable to the relative change in the parameter (Chitnis et al., 2008). The normalized sensitivity index of R 0 , that depends differentiably on a parameter p, is defined as The normalized sensitivity indices of R 0 are summarized in Figure 3. For an example of the calculation, see Appendix A.2. Results show that the biting rate of mosquitoes, b, is the most impactful parameter in the computation of R 0 . That is, an increase in the number of times per day a single mosquito bites humans will correspond to the greatest increase in the probability of humans becoming infected with Zika. The next most sensitive parameters are μ −1 V and γ −1 , the lifespan of vectors and the infectious period for humans, respectively. Increasing μ −1 V results in a longer lifespan of mosquitoes, thus increasing the time mosquitoes are able to transmit the disease and as a consequence increasing R 0 . Increasing γ −1 results in a slower recovery rate, and thus infectious humans are able to transmit the disease for a longer period of time, and as a result R 0 increases.
The least sensitive parameter is q, the proportion of humans who are symptomatic. This is not surprising as sexual transmission makes up only 14% of R 0 . However, if we restrict our analysis to the effect that asymptomatic individuals can have on R 0S , the effect of the behaviour of asymptomatic individuals is significant, see Section 4 and Figure 4.

Wolbachia as a form of vector control
The Model presented here, as many mathematical epidemiological models, is used to experiment with control measures that can be used by governments to decrease the size and length of a Zika epidemic. In this work, the Wolbachia bacterium is introduced into the vector population as a measure of control. As discussed in Section 1, Wolbachia bacterium can play a major role in reducing the mosquito population and preventing the spread of Zika. Previous mathematical models that include Wolbachia do so by explicitly modelling Wolbachia infected mosquitoes (Hughes & Britton, 2013;Ndii et al., 2012). The model presented here takes a simpler approach.
In order to introduce the effect of control, the model focuses on the parameters affected by Wolbachia. When male mosquitoes are infected with Wolbachia, sperm and eggs are unable to form viable offspring, thus lowering the value of λ V . While this is important, the parameter λ V does not appear in the expression of R 0 . To calculate the control measures, a scaling parameter is attached to λ V and is then used to determine the time course of E H and E V with a 25%, 50% and a 75% reduction in vector births (see Figure 5). The dark grey dot-dashed curve is a 25% reduction. The grey long-dashed curves are 50% and the light grey-dashed curves are a 75% reduction: (a) exposed humans and (b) exposed vectors.
Another way Wolbachia can be used as a control measure is in a reduction in transmission from vectors to humans. A study was done in 2010 that found at least 37.5% of female mosquitoes that had been infected with both Dengue and Wolbachia were unable to transmit the disease (Bian, Xu, Lu, Xie, & Xi, 2010). If a similar effect is assumed with Zika, the transmission probability, a H , from an infectious mosquito to a susceptible human is decreased. The only entry in the NGM K that is targeted with the control strategy of decreasing a H is the element K 12 . Thus we calculate the target reproduction number T S where S = {(1, 2)} as described in van den Driessche (2017). For information regarding the type reproductive number, calculated when the control strategy is aimed at particular host types only, see Heesterbeek and Roberts (2007); Roberts and Heesterbeek (2003). For more detailed work on the target reproduction number, calculated when control can be targeted at interactions between types, see Shuai, Heesterbeek, and van Den Driessche (2013).
Thus, following the recipe in van den Driessche (2017), the target matrix, K S , is Note that ρ(K − K S ) = R 0S < 1 and thus the target reproductive number T S is given by If the transmission from vectors to humans can be reduced by a fraction of at least 1 − 1/T S = 1 − (1 − R 0S )/(R 0V ) 2 = 0.6871 then Zika will die out. We also provide a similar figure as for λ V by attaching to the parameter a H a scaling parameter and seeing again what a 25%, 50% and 75% reduction in transmission from vectors to hosts does to the system (see Figure 6). Control measure for reducing vector transmission to humans. The solid black curves represent no reduction. The dark grey dot-dashed curve is a 25% reduction. The grey long-dashed curves are 50% and the light grey dashed curves are a 75% reduction: (a) exposed humans and (b) exposed vectors.

Results
The model presented is used to determine the effect of the sexual behaviour of asymptomatic individuals in the spread of Zika. Figure 4 shows the percentage change in R 0S , the basic reproductive number tied to the human-to-human disease transmission, as ψ A , product of the sexual infectiousness and the rate at which asymptomatic individuals engaged in unprotected sex, is varied. As the amount of unprotected sexual encounters had by asymptomatic individuals increases, so does the probability that Zika is spread, which should not be surprising. What is more telling is the linear relationship between ψ A and R 0S , and how doubling ψ A results in almost doubling R 0S . The model was also used to test Wolbachia as a vector control strategy. Using the parameter values from Table 2, scaling parameters are used to give a 25%, 50%, and 75% reduction in the parameters λ V and a H . The model is run until equilibrium and the time course of the Exposed classes are plotted to show the results. Figure 5 displays the results for scaling the vector birth rate, λ V . The plots show that as the vector births are reduced, the length of the outbreak gets shorter as does the magnitude of the outbreak. For the vector transmission to humans, Figure 6 shows the magnitude of the outbreak is reduced, however, the length of it is prolonged. The outbreak is not as severe but lasts for an extended period of time. The length is much longer for the 50% reduction, but at the 75% reduction the outbreak dies out very quickly. Furthermore, the target reproduction number that was found in Section 3 showed that if a H is reduced by a fraction of at least 0.6871, then R 0 will be less than one. This can be seen in Figure 7 where R 0 graph crosses the horizontal line y = 1 when a H reaches the value of 0.3 ≈ (0.6871)(0.4240).

Discussion and future work
The model presented here aims to simulate an outbreak of the vector borne Zika virus taking into consideration symptomatic and asymptomatic humans and the possibility of sexual transmission. In order to gain better understanding of the spread of the disease, the basic reproductive number R 0 was computed and separated into the two components, Figure 7. Graph representing R 0 as a function of the vector to human transmission rate. In this graph, it can be observed that if the transmission from vectors to human can be reduced by a fraction at least 0.6871, then the vector-host disease will die out. This information can be used to evaluate any vector control strategy, in particular, the use of Wolbachia.
the basic reproductive number due to vector to human transmissions and the basic reproductive number due to human to human transmissions, R V and R 0S respectively.
Given that 80% of individuals do not present any symptoms, the model was used to explore the effect the sexual behaviour of asymptomatic individuals can have on the spread of the disease by plotting the change in R 0S against changes in the ψ A , the rate of sexual transmission of asymptomatic individuals and showing the linear relationship between these two variables, R 0S and ψ A .
The model was then used to investigate the effect of implementing vector control strategies, in particular the use of Wolbachia focused on the two types of control: the reduction in vector births due to cytoplasmic incompatibility where the embryos do not develop; and through the reduction in vector transmission to humans. The focus on vector transmission seems to not be as effective as it requires an almost 70% reduction in order for R 0 to be reduced below 1. It also seems to increase the length of the outbreak considerably. The vector birth control seems to be more effective. The model shows that even a 50% reduction in vector births can shorten the outbreak time and reduce its magnitude.
It is important for the projects that are releasing Wolbachia infected mosquitoes to monitor the reduction in vector births. Some projects are only releasing male Wolbachiainfected mosquitoes. However, females are less likely to mate with Wolbachia-infected mosquitoes and because of cytoplasmic incompatibility, Wolbachia will not continue to strain to future broods (Jiggins, 2017). Wolbachia-infected females have a selective advantage over uninfected females because they have a normal brood size, regardless of the Wolbachia-infection status of the males they mate with Sullivan and O'Neill (2017). This advantage would help Wolbachia spread through (and ultimately reduce) the mosquito population: the male offspring from a Wolbachia-infected female would essentially be sterile, and the female offspring would continue the Wolbachia strain to future broods.
There are several directions in which the work presented here can be expanded. A simple one is to use the model presented here to test other measures of control. Similar models have tested different controls measures such as the use of bed nets to reduce vector-host infections or the use of condoms to decrease sexual transmission among others. Another direction that can be taken is to split the human population by gender and consider the different recovery rates between females and males and all different sexual transmission possibilities: males to males, male to female, female to male and female to female further complicating the model. This direction seems particularly important to take as the recovery rate for human females is on the order of days while Zika has been shown to remain in male semen up to six months after infection. We believe these explorations will lead to interesting and innovative dynamics as Zika seems to be the only vector born disease to exhibit such behaviour.
We let (x 1 , x 2 , x 3 , x 4 , x 5 ) = (E H , I HS , I HA , E V , I V ) and write the Jacobian matrix of the infection subsystem as T + , where T is the transmission matrix and the transition matrix. Hence T contains the entries where an epidemiological birth occurs, and all other epidemioligcal events are incorporated in the model via . Thus Referring to infected states with indices i and j, where i, j = 1, . . . , 5, the ijth entry of T is the rate at which individuals in infected state j have epidemiological offspring in infected state i for the linearized system. Thus we note that rows two, three, and five, of T consists of zeroes only as individuals in infected state j = 1, . . . , 5 do not produce new cases in states I HS (row 2), I HA (row 3) and I V (row 5). The transition matrix , corresponding to all other changes of state, is given by Thus the NGM with large domain, K L = −T −1 , is We note that although −1 is not difficult to compute, it is also easily found by using the biological interpretation of − −1 . That is, the ijth element of − −1 is the expected time that an individual presently in infected state j will spend in state i as the disease runs its course. Thus the diagonal terms of − −1 are the number of days a human or vector spends in its present state. As a more interesting example, a vector in infected state E V is expected to spend β V /(β V + μ V ) × 1/μ V days in state I V , where the first factor is the probability that an individual vector survives state E H and moves to state I V , and the second factor is the expected amount of days it will spend in state I V .
Given the epidemiological interpretations of T and − −1 , the ijth entry of K L is the expected number of state-i-infected offspring an individual presently in infected state j will produce throughout its future infected life. Because the infected states I HS , I HA and I V are not states-at-infection, the matrix K L has rows two, three and five, exactly zero. Thus to obtain K from K L we let the auxiliary matrix, E, be given by Then, the NGM is