Transmission dynamics of Zika virus with multiple infection routes and a case study in Brazil

The Zika virus (ZIKV) is a serious global public health crisis. A major control challenge is its multiple transmission modes. This paper aims to simulate the transmission patterns of ZIKV using a dynamic process-based epidemiological model written in ordinary differential equations, which incorporates the human-to-mosquito infection by bites and sewage, mosquito-to-human infection by bites, and human-to-human infection by sex. Mathematical analyses are carried out to calculate the basic reproduction number and backward bifurcation, and prove the existence and stability of the equilibria. The model is validated with infection data by applying it to the 2015–2016 ZIKV epidemic in Brazil. The results indicate that the reproduction number is estimated to be 2.13, in which the contributions by mosquito bite, sex and sewage account for 85.7%, 3.5% and 10.8%, respectively. This number and the morbidity rate are most sensitive to parameters related to mosquito ecology, rather than asymptomatic or human-to-human transmission. Multiple transmission routes and suitable temperature exacerbate ZIKV infection in Brazil, and the vast majority of human infection cases were prevented by the intervention implemented. These findings may provide new insights to improve the risk assessment of ZIKV infection.


Modeling framework
A mechanistic framework for simulating ZIKV transmission patterns is established in this section.Inspired by recently-developed mathematical models 9,[16][17][18][19] , and based on epidemiological feature of ZIKV, a new dynamic system is proposed, which takes into account human-to-human (sexual transmission) and human-to-mosquito (bite transmission and sewage transmission) interactions by using compartmental and deterministic principle.The total numbers of human, larval and adult mosquitoes are represented by N h , M v and N v .Given the stability of human demographics and the extremely low rate of ZIKV-related mortality, N h is assumed to be a constant.
Based on the characteristics of ZIKV infection, it is further assumed that: (1) the humans are divided into five categories: susceptible S h , latent E h , symptomatic infected I s , asymptomatic infected I a and recovered R h ; (2) the larval stage can be divided into two categories: uninfected A v (including eggs, larvae and pupae) and infected J v ; (3) the adult female mosquitoes are divided into three categories: susceptible S v , latent E v and infected I v .
It is known that mosquito oviposition is linked to mature females and the ability of mosquitoes to develop oviposition habitats.If there are too many eggs in the oviposition habitat or too few nutrients and water resources, the females will lay fewer eggs or choose another site 21,22 .In addition, larvae and pupae need water or nutrients complete their development 21 .This biological phenomenon can be expressed by a Logistic model which explicitly incorporates the idea of limited carrying capacity resources 21,22 .Hence, the per capita oviposition rate is given by θ(1 When ZIKV invades an area, humans and mosquitoes there could be infected with certain probability.It is governed by the following rules.A susceptible human may be infected with ZIKV from the bite of infectious mosquitoes at rate b or through human contact with infected people at rate β .Infected humans undergo an incu- bation period 1/δ .After that they could be symptomatic φ and asymptomatic 1 − φ , and then the former follows by an infectious period of mean duration 1/η until recovery.The larvae may be directly infected by the virus sewage at rate p and grow up to become infectious adult mosquitoes at rate α .Susceptible mosquitoes become infected at rate a by biting infectious humans, and the infected mosquito experience an extrinsic incubation period 1/γ , that is followed by an infectious state from which they do not recover.The mean larvae lifespan and adult lifespan are 1/f and 1/µ , respectively.Flow chat is shown in Fig. 1.
Accordingly, the following ODEs are used to simulate the transmission dynamics of ZIKV: Detailed explanations of the variable and parameters are shown in Table 1.All the variables and parameters are non-negative.Similar to the proofs in 21,22 , it is obtained that the following is positively invariant and attracting under the flow described by the system (1).
Temperature, as the most important factor in the modulation of mosquito-borne diseases, is considered by including in the model parameters, where the parameters are functions of temperature T. Based on observations from laboratory studies [23][24][25] , they are presented below. (1)

Mathematical analysis
The transmission dynamics of the proposed model are analyzed mathematically, in which the epidemic threshold (i.e., the basic reproduction number) is calculated by the next generation matrix method, and the evolution trends of the model solutions are investigated by stability theory.

Basic reproduction number
The basic reproduction number R 0 is interpreted as the average number of secondary cases that are produced by a single primary case in a fully susceptible population, acting as the critical measure of the transmissibility 31 .The basic reproduction number R 0 is calculated by using the theory of next generation matrix 31 .It is written as 31 : R 0 = ρ FV −1 , where F is the rate of occurring new infections, and V is the rate of transferring individuals outside the original group.Here ρ represents the spectral radius of matrix.
Let the right hand side of the system (1) equal to zero, it is obtained the disease-free equilibrium of the model where It follows from the next generation matrix method 31 that Using R 0 = ρ FV −1 31 , the basic reproduction number is calculated as follows: (6) a(T) = max 0.001044T(T − 12.286)(32.461− T) 1/2 , 0 ,  www.nature.com/scientificreports/where Some typical features can be observed from Eq. ( 9): (a) R hh 1 and R hv 1 represent the parts of the basic reproductive number R 0 contributed by sexual transmission and mosquito-borne transmission, respectively; (b) R hh 1 is determined by human-to-human transmission; (c) R hv 1 is determine by parameters related to mosquitoes and humans, which can be further divided into two stages of larval transmission and adult transmission.

Model stability
Mosquitoes breeding in contaminated water can become infected with ZIKV, which can shorten the virus transmission cycle and promote the spread of ZIKV 14 .Based on this, the following mathematical analysis is divided into two cases, with or without larval infection in the virus sewage, corresponding to p = 0 or p = 0. p = 0 (without larval mosquito infection in the virus sewage) In this case there is no larval infection J v item in the system (1).The basic reproduction number of the available model is where Theorem 3.1 When p = 0 , if R 0 < 1 , then the disease-free equilibrium of the system (1) is locally asymptotically stable.
The proof is similar to that of Theorem 3.5, so it is omitted here.The expression of endemic equilibrium is denoted by Letting the right-hand side of system (1) to be zeros, direct calculation yields an equation about infectivity * h as Let

It is further obtained that
where It follows from Eq. ( 11) that where Here * h can be obtained by solving the quadratic expression (12), which also determines the expression of the endemic equilibrium.When .
Vol:.( 1234567890 www.nature.com/scientificreports/ In this case, Since c 0 > 0 and c 2 < 0 (when R 0 > 1 ), it is obtained a unique solutions for Eq. ( 12), which indicates that the system (1) has a unique endemic equilibrium when R 0 > 1 .Knowing that c 2 > 0 when _ R 0 < 1 , there are three scenarios under this condition: h + c 2 is on the negative half axis of x, and there is no intersection point between f * h and the positive half axis of x, so the system (1) has no endemic equilibrium.(ii) If c 1 < 0 and c 2 1 = 4c 0 c 2 , there is only one intersection point between f * h and the positive half axis of x, so the system (1) has a positive equilibrium point.(iii) If c 1 < 0 and c 2 1 > 4c 0 c 2 , f * h has two intersection points with the positive half axis of x, then the system (1) has two positive equilibrium points.
The above analysis concludes the following theorem.Theorem 3.2 When p = 0 , if R 0 > 1 , then c 2 < 0 , the system (1) has a unique endemic equilibrium.If R 0 < 1 , the following conclusions can be drawn: then System (1) has two endemic equilibria.
µ , then the system (1) has a backward bifurca- tion; if R 0 > 1 , the system (1) has a forward bifurcation and the endemic equilibrium is locally asymptotically stable.
T to be the left-hand and righthand sides of system (1).Thus system (1) can be written as The Jacobian matrix of system (1) at the disease-free equilibrium E 0 is The probability of mosquito-to-human transmission b is selected as the branching parameter.When R 0 = 1 , the critical value b is solved, Therefore, when b = b * , the characteristic equation of the system is w h e r e h( ) Vol 2 +a 2 a 3 +a 2 a 4 +a 3 a 4 −H 3 ) > 0 .Hence, it follows from Rrouth − Hurwitz theorem that the real parts of the characteris- tic roots of the equation h( ) = 0 are all negative.
Furthermore, by calculating the eigenvalue of the matrix J b * (E 0 ) , it is found that J b * (E 0 ) has a zero eigenvalue, and other eigenvalues have negative real parts.Therefore, the system (13) satisfies the conditions of the Central Manifold theorem.The right eigenvector and the left eigenvector corresponding to the eigenvalues = 0 of the matrix J b * (E 0 ) are respectively denoted as where w 1 = −a 4 w 2 /d , w 2 > 0 , w 3 = φδw 2 /a 3 , w 4 = (1 − φ)δw 2 /a 3 , w 5 = ηδw 2 /(da 3 ) , w 6 = 0 , w 7 = a 2 w 8 /µ , w 8 = µw 9 /γ , w 9 = w 2 (a 4 − H 3 /a 3 )/(bc) , and According to Castillo-Chavez and Song theorem 32 , it is calculated The sign of B will determine the occurrence of a backward bifurcation in a given model.When b = b * , the positive and negative of the coefficient A determines the local dynamical properties of the disease-free equilibrium 32 .Hence, it follows from Castillo-Chavez and Song theorem that system (1) undergoes a backward bifurcation at 3 a 2 4 µ .When R 0 > 1 and close to 1, it has A < 0 and thus the system has a forward branch and the endemic equilibrium is locally asymptotically stable.Theorem 3.4 When p = 0 , if R 0 > 1, the endemic equilibrium of system (1) is globally asymptotically stable.

Proof Define the following functions
Using the inequality 1 − x + ln x ≤ 0 , for x > 0 , differentiation yields: Vol:.( 1234567890 www.nature.com/scientificreports/With the constants a ij above and A = [a ij ] , we construct the (strongly connected) directed graph Ŵ(A) in Fig. 2. If and only if a ij > 0 , there is a weighted arc i, j .Along each of the cycles on the graph, one can verify that G ij = 0 , for instance, G 5,4 +G 1,5 +G 2,1 +G 4,2 =0 , and so on.Then, by Theorem 3.5 in 33 , there exist constants c i such that Q = i c i Q i is a Lyapunov function for system (1).To find the constants c i , we use the combinatorial Then the function For the left subsystem, www.nature.com/scientificreports/One can show that system (19) has a unique equilibrium R * h , A * v , and that this point is globally asymptotically stable for this system.Therefore, the largest and only invariant set in S is the endemic equilibrium Using LaSalle¡¯s Invariance Principle, we conclude that the endemic equilibrium E * 0 globally asymptotically stable in .
p = 0 (mosquitoes hatched in contaminated water can be infected with ZIKV) Theorem 3.5 When R 0 < 1 , the disease-free equilibrium of the system (1) is globally asymptotically stable.

Let
It can be seen that the right side of the system (20) is the right side of the matrix F − V .According to R 0 = ρ FV −1 < 1 , it can be seen that the system (20) has only a balance point E h , I s , I a , J v , E v , I v =(0, 0, 0, 0, 0, 0) .Therefore, every non-negative solution of (20) satisfies Because the system (20) is linear, the disease-free equilibrium of the system (20) is globally asymptotically stable.
According to the comparison theorem

Thus
The disease-free equilibrium is globally attractive, and it is locally stable.So the disease-free equilibrium of the system (1) is globally asymptotically stable.Theorem 3.6 When R 0 > 1 , the system (1) exists endemic equilibrium.Substituting * h by 0 and +∞ , it follows that f (+∞) = −1 < 0 , and According to the existence theorem of zero point, when R 0 > 1 , there exists * h > 0 , such that f * h = 0 , indicating the existence of endemic equilibrium for the system (1).

Model application
The massive outbreak of Zika in Brazil during 2015 and 2016 is used as a typical prototype to validate the proposed model.The estimate_R function in EpiEstim package of R language is used to calculate the effective reproduction number R t based on weekly morbidity time series, intergenerational distribution and window period 34 .R t is the average number of people someone infected at time t can infect over their entire infectious lifespan, which can quantify the immediate transmissibility.Here R t is used to determine the parameter b in the model.Since clinical studies have shown that the viral load of asymptomatically infected patients with Zika is about half that of symptomatic patients, it is assumed that the transmission probability of asymptomatic people is half that of symptomatic people 26,27 .In addition, due to the large temperature difference between the north and south of Brazil, the average weekly temperature T = 24.15• C of the three cities of Manaus, Brasilia and Porto Alegre (divided into the north, central and south of Brazil) is employed in the model parameters expressions.
In order to verify the model of case study, MCMC method is used to quantify the uncertain parameters: ϑ , β and p.In the absence of surveillance data, it is assumed that the initial values of total mosquitoes sizes are at the positive stable level, M v (0) and N v (0) are given by ().The initial values of the numbers of vector in latent and infected states are also estimated by MCMC method.The MCMC is run for 200,000 iterations for each parameter and the posterior distributions are compiled from the final 80% of the iterations.The model is validated with 95% confidence interval of posterior estimations by sampling 1000 times of the posterior distributions of the estimated parameters and by incorporating them into the model.The normalized forward sensitivity and global sensitivity are used to quantify the importance of parameters related to the modeling incidence and the reproduction number R 0 .The normalized forward sensitivity is computed by ∂R 0 ∂x x R 0 for parameter x.The global sensitivity is realized by computing partial rank correlation coefficient (PRCC) based on Latin hypercube sampling for the model inputs and outputs 35 . Vol.:(0123456789)

Result
The reporting data shows that Brazil witnessed the first human infection in January 2015, but few case is recorded in the following 2 months.ZIKV began to expand since April, with the first epidemic peak in early July.After a continuous low incidence, human infections reached a higher peak in February 2016 and then decrease rapidly until zero in September.By fitting the epidemic curve of the cumulative number of weekly cases by the model, the estimated results of unknown parameters are shown in Table S1 and Fig. S1.As shown in Fig. 3a, it is observed that the estimated parameters enable the model to draw a good fitting capacity, except the first peak.The fitting deviation is possibly due to the temporal heterogeneity of transmission parameters and detection efficiency of human cases.Uncertainty analysis indicates that the model is robust in exploring transmission dynamics, which can draw consistent evolution of weekly accumulate incidence in case of random sampling (see Fig. S1 in SI).
Sensitivity analysis is used to quantify the response of model outputs to parameter variation.Results from Fig. 3 (b) indicate that the most sensitive parameter to the modeling infection is environmental capacity rate of mosquitoes ( ϕ ), followed by the human-to-human transmission rate ( β ) and the initial value of infected adult mosquitoes ( I v ).Yet the model output is not sensitive to the transmission rate from infected person fecal to larval (p).
The contributions of different routes to the ZIKV infection in Brazil is estimated by splitting the values of the basic reproduction number R 0 .Substituting the estimated parameters (in which b is chosen as the average value of the previous 7 weeks) into Eq.( 9) yields the basic reproduction number to be R 0 = 2.13 (95% CI 1.61-2.64), in which the contribution of mosquito transmission is 2.04 (95% CI 1.53-2.55).Moreover, the transmission of the virus by mosquitoes could be divided into adult mosquitoes and larval mosquitoes in sewage, in which the latter contribution to R 0 is estimated to be 0.48 (95% CI 0.26-0.71).
Results of sensitivity analysis toward R 0 are presented in Fig. 4. Two different methods of sensitivity analysis show similar results, indicating the robustness of R 0 to the parameters used.It is observed that the mortality of adult mosquitoes ( µ ), infection period of human beings ( η ), mosquito biting rate (c), environmental capacity rate of mosquitoes ( ϕ ), human-mosquito transmission rates (a and b) and transition rate of mosquitoes from larval to adult ( α ) are most sensitive parameters to determine R 0 .Yet R 0 is less sensitive to the transmission rates from fecal to larval and from human to human (p and β ), oviposition rate of adult mosquitoes ( θ ) proportion of symptomatical infections ( φ ), human death rate (d) and relative infectivity of asymptomatic infections (q and τ ).These parameters could play minor roles in causing human infection of ZIKV.
Figure 5 shows the effects of different transmission paths of the Zika outbreak in Brazil on the cumulative number of human infections from January 2015 to September 2016.The fitting results demonstrate that the total infections could be 304,648 (95% CI 304,096-305,199).If without human infection by sex or without larvae mosquito infection by sewage, this number could be 242,487 (95% CI 242,025-242,949) or 297,115 (95% CI 296,558-297,674).If without both of the above transmission routs, it could become 236,455 (95% CI 235,994-236,918).Of these three assumptions, the total number of human infections decreased by 20.4%, 2.47% and 22.38%, respectively.Moreover, in the absence of asymptomatic infection, the cumulative number of human infections could be 254,715 (95% CI 254,238-255,192), making human infections decrease by 16.39%.In this case, the above-mention three circumstance lead to the declines of human infections by 6.43%, 0.79% and 7.17%, respectively.Figure 6 shows the time evolution of human infection in cases of different human-to-human transmission rate (β) , larval mosquito transmission rate in sewage (p) and temperature.It is observed that the increase of transmission rates through sex and sewage would lead to moderately higher number of human cases and slightly faster of disease transmission, Such effect caused by high sex transmission rate is more significant.Furthermore, temperature around 28.7 • C is most favourable for ZIKV infection, and temperature away from this value could cause low morbidity and low infection risk.
Figure 7 shows the effectiveness of the implementation of control measures on the spread of ZIKV infection in Brazil during January 2015 and November 2016.The intervention is measured by the effective reproduction number R t .Here the values of R t are fixed to be 2.17, 2.12, 1.98, 1.83, 1.70 and 1.62, which are its average values of in the previous 1-5, 1-10, 1-15, 1-20, 1-25 and 1-30 weeks, respectively.The simulation results indicate that (1) if R t = 2.17 (few intervention measures), the cumulative number of human infections could reach 37.2 mil- lion, that is over 121 times of reported cases, with early peak of new cases as 3.7 million around the 37th week; (2) if R t = 1.62 (limited intervention measures), the cumulative number of human cases could be 36 million, with peak number of new cases as 2.1 million around the 57th week.It is observed that the decrease of R t causes an evident decline of humans infections and quick arriving of peak infection.

Discussion
A modeling framework for inferring ZIKV transmission patterns is attempted in this paper.Technologically, a new SEIR-AJ-SEI dynamic model is established, which couples the ZIKV circulation among/between mosquitoes and humans under potential routes, including mosquito bite, sex contact, and sewage breed.Compared with existing ZIKV model, such as using a discrete stochastic SEIR-SEI model to predict the optimal effect of bednets, infection treatment and insecticide spraying on disease transmission 18 , using a SEIIAR-SEI model to infer the impact of mosquito-borne and sexual transmission on ZIKV spread 9 , using a continuing climate-driven SEIIIR-SEI model to study threshold dynamics in a seasonal model of Zika virus disease 19 , using a high-dimensional ODE system to describe the joint dynamics of Zika and dengue 17 , and using partial differential equations model to understand how spatial heterogeneity of the vector and host populations influence ZIKV dynamics 16 , the proposed model is a combined update of recent works, which converges more dynamical detail.
First, the ZIKV transmission dynamics are clarified mathematically.The basic reproduction number R 0 is calculated by using the next generation matrix, which is found to be determined by all the transmission routes.It is verified that the disease-free equilibrium is locally stable when the associated basic reproduction number is less than unity, and there exits endemic equilibrium when basic reproduction number is larger than unity.If without larval infection in virus sewage, central manifold theorem demonstrates that the system is capable of undergoing the phenomenon of backward bifurcation, under which the stable disease-free equilibrium would co-exist with a stable endemic equilibrium and an unstable endemic equilibrium.Such phenomenon could be caused by the combined effects of multiple transmission routes.The existence of larval infection in virus sewage would cause more complicated transmission dynamics.Dynamic details indicate that great efforts could be needed for preventing ZIKV infection.
Second, the proposed model is further validated to explore more transmission dynamics by fitting the reported cases of ZIKV infection in Brazil from January 2015 to September 2016.Two important insights based on this study may provide scientific clues for evaluating the infection risk and guiding control.
On the one hand, the impact of the transmission routes on the Zika epidemic is evaluated.It is found that mosquito-human infection by bite is still the prominent path for ZIKV occurrence, accounting for 85.7% of the basic reproduction number R 0 , but the contributions by sexual transmission and larval transmission in sewage are 3.5% and 10.8%, respectively, both of which are less than 1.Therefore, these two routes of infection are not enough to trigger a large-scale outbreak of ZIKV, which is consistent with previous studies 9,11,15 .The sensitivity analysis further confirms that ZIKV infection is dominated by the parameters related to mosquito ecology, rather than those parameters related to asymptomatic or human-human transmission.Yet the latter two transmission routes can limitedly accelerate the development of Zika epidemic and prolongs the outbreak time.Such effects   www.nature.com/scientificreports/would be more significant for higher concentration of ZIKV in sewage and larger probability of human-tohuman infection.Therefore, the prevention and control of ZIKV should target at reducing the infection through mosquito bite, and do not ignore the infection through sex and sewage.On the other hand, the situations of ZIKV transmission in Brazil are evaluated.The present study suggests that multiple modes of transmission and suitable temperature may be responsible for the large outbreak of ZIKV in Brazil in 2015-2016.Nevertheless, the intervention implemented in Brazil plays an important role in controlling ZIKV infections.If without intervention, the number of human infections with ZIKV in Brazil would increase rapidly 37 , and result in more than 37.2 million cases.After the implementation of control measures in May 2015, the number of effective reproduction R t decreased from 1.69 to less than 1, leading to a rapid decline in morbidity.However, human infections began to rebound rapidly in November of the same year, leading to the second peak of cases in Brazil at 59th week.Meanwhile, control measures were intensified, including mobilizing Brazilian army to support community health services 37 , house-to-house visits and the elimination of potential Aedes breeding sites 37 , aerial spraying of the product to kill larvae or adult mosquito 38 , and reduction of breeding sites through drainage of standing water, waste management and education about mosquitoes and personal protection measures 39 .Under these control measures, the incidence and effective reproduction number dropped fast.The present study indicates that such comprehensive and intensified ZIKV control strategies are highly effective in curtailing ZIKV transmission.
The following limitations need to be clarified.First, since the proposed model is only used to fit the surveillance data reported in Brazil, there may be geographical differences in its application to other countries.Second, since it is impossible to obtain the true values of some model parameters, they are extracted from the literature or are estimated by MCMC method.Third, the study dose not include all the underlying factors, such as human mobility and climate (except temperature).However, the model takes into account the most influential factors and incorporates model parameterizations, providing confidence in the model output for future analysis.

Figure 1 .
Figure 1.Flow diagram of the ZIKV transmission among humans and mosquitoes.The black solid lines indicate the progress of infection and ecology.The yellow lines indicated the infection pathes.

Figure 3 .
Figure 3. (a) The fitting results of the ZIKV cases in Brazil from the proposed model.The light colored area is the 95% confidence intervals (CIs) for all 1000 simulations.(b) Global sensitivity analysis for weekly incidence, where the PRCCs are the mean values in each week, and * indicates a significant difference from zero (with p value < 0.01).

Figure 4 .Figure 5 .
Figure 4. Sensitivity analysis of the basic reproduction number R 0 .(a) Normalize forward sensitivity; (b) global sensitivity with PRCC.

Figure 6 .
Figure 6.Impacts of (a) sexual transmission rate β , (b) sewage transmission rate p and (c) temperature on the number of new infections cases, which is achieved by simulating the proposed model with other parameters equal to those of the fitting results.

Figure 7 .
Figure 7. Impacts of different effective reproduction number on the number of human infection in Brazil during January 2015 and November 2016, which is achieved by simulating the proposed model with other parameters equal to those of fitting results.

Table 1 .
Description of parameters in the model (the time unit is week or per week).[a] b is determined by the effective reproduction number R t in the case study.[b] β and p are estimated using the MCMC approach in the case study.[c] K is assumed to be proportional to population size, i.e.,K = ϕN h , where ϕ is estimated by MCMC method.
36azil is located in the eastern South America, with longitudes 35 • W to 74 • W and latitudes 5 • N to 35 • S. It is the largest country in South America, whose area is about 8,514,900 square kilometers.Brazil has typical tropical climate, with the annual average temperature as 20-28 • C. Such climate is suitable for the growth of Aedes mosquitoes.As a result, Brazil faces public health problems related to mosquito-borne diseases.The study use medical records of human Zika infections reported in Brazil from January 2015 to September 2016.The weekly numbers of cases is collected from the World Health Organization36and the Brazilian Ministry of Health, which is shown in TableS2(see Supplementary Information [SI]).Brazilian population data is obtained from the Brazilian Institute of Statistics (https:// www.ibge.gov.br/).Brazil's weekly average temperature record is extracted from the Brazilian Meteorological Agency (https:// previ sao.inmet.gov.br/).