Multiple instances of Zika virus epidemic have been reported around the world in the last two decades, turning the related illness into an international concern. In this context the use of mathematical models for epidemics is of great im- portance, since they are useful tools to study the underlying outbreak numbers and allow one to test the e ectiveness of di erent strategies used to combat the associated diseases. This work deals with the development and calibration of an epidemic model to describe the 2016 outbreak of Zika virus in Brazil. A system of 8 di erential equations with 8 parameters is employed to model the evolu- tion of the infection through two populations. Nominal values for the model parameters are estimated from the literature. An inverse problem is formulated and solved by comparing the system response to real data from the outbreak. The calibrated results presents realistic parameters and returns reasonable de- scriptions, with the curve shape similar to the outbreak evolution and peak value close to the highest number of infected people during 2016. Considera- tions about the lack of data for some initial conditions are also made through an analysis over the response behavior according to their change in value. Keywords: Zika virus dynamics; nonlinear dynamics; mathematical biology; SEIR epidemic model; model calibration; inverse problem Corresponding author. Email addresses: eber.paiva@uerj.br (Eber Dantas), michel.tosin@uerj.br (Michel Tosin), americo@ime.uerj.br (Americo Cunha Jr) Preprint submitted to Applied Mathematics and Computation July 4, 2018 1. Introduction The Zika virus is a avivirus that upon infection in humans causes an illness, known as Zika fever, identi ed commonly with macular or papular rash, mild fever and arthritis [1, 2]. It is mainly a vector-borne disease carried by the genus Aedes of mosquitoes [2, 3], while in a lesser amount it is also transmitted via sexual interaction [4, 5], and contamination by blood transfusion is under investigation [6]. The Zika virus was rst isolated in primates from the Zika forest in Uganda in 1947 [7]. Evidences of the virus in humans were found in Nigeria in 1968 [8]. An epidemic occurred in 2007 on Micronesia [9], followed by multiple outbreaks on several Paci c Islands between 2013 and 2014 [10, 11]. The rst Zika virus autochthonous case in Brazil was reported around April, 2015 [12], and nearly 30 ;000 cases of infection were already noti ed by January 30, 2016 [13], along with the Pan American Health Organization being informed in the same month about locally-transmitted cases on numerous continental and island territories of America [14]. The Brazilian Ministry of Health registered 215;319 probable cases of Zika fever (130 ;701 of which were con rmed) until the 52th epidemiological week (EW) of 2016 [12]. The Zika epidemic has been causing concern in the international medical community, health authorities and population, specially due to an association between the Zika virus and other diseases, such as newborn microcephaly [4, 15] and Guillain-Barr e syndrome [16], whose correlation to the Zika virus was considered by the World Health Organization a \scienti c consensus" [17]. In this epidemic scenario, the development of control and prevention strate- gies for the mass infection is a critical issue. A mathematical model capable of providing a description of the infected people throughout an outbreak is an useful tool that can be employed to identify e ective and vulnerable aspects on disease control programs [18, 19, 20]. Furthermore, for an epidemic model to be truly useful it must undergo a judicious process of validation [21, 22], which con- sists in comparing model predictions with real data in order to evaluate if they are realistic. In general, the rst predictions of a model do not agree with the observations, possibly due to inadequacies in the model hypotheses or because of a poor choice for the model nominal parameters. The rst case invalidates the model, but the second can be amended through a procedure known as model calibration , where a set of parameters that promote a good agreement between predictions and observations is sought. This work is one of the results in a rigorous ongoing process of identi cation and validation of representative models to describe Zika virus outbreaks in Brazil [23, 24]. For this purpose, a SEIR-SEI mathematical model is adapted to the Brazilian scenario. This speci c SEIR-SEI description has been successfully used before for the outbreaks in Micronesia [25] and French Polynesia [26]. Some assumptions were also based on similar studies performed over SEIR dynamical systems [27, 28, 29, 30]. The nominal values of the model parameters belong to characteristics of the Zika infection and its vector, quantitatively estimated in the literature or published by health organizations. Predictions are obtained from numerical simulation and further heuristic manipulation, followed by a 2 comparison to real data of the outbreak as an initial e ort to validate the model. In sequence, a rigorous process of model calibration is performed through the formulation and solution of an inverse problem. The rest of this paper is organized as follows. In section 2, the mathematical model is described and estimation of nominal values for the model parameters is discussed. In section 3, the forward and inverse problems are formulated and solved, where results are detailed and a subsequent comparison between model predictions and experimental data is made to calibrate the model. Finally, in section 4, the main contributions of this work are emphasized and a path for future works is suggested. 2. Epidemic model for Zika virus dynamics 2.1. Model hypotheses This work utilizes a variant of the Ross-Macdonald model [31] for epidemic predictions, separating the populations into a SEIR-SEI framework (susceptible, exposed, infectious, recovered) [32, 33, 34]. Each category represents the health condition of an individual inside such group at time t, with respect to the con- sidered infection. The susceptible group, denoted by S(t), represents those who are uncontaminated and are able to become infected. The exposed portion of the population, E(t), comprehends anyone that is carrying the pathogen but is still incapable of transmitting the disease. While the infectious individuals, I(t), can spread the pathogen and may display symptoms associated with the illness. Finally, the recovered group, R(t), contains those who are no longer infected. The populations of humans and vectors are segmented into the aforesaid classes (excepting the group of recovered vectors), as Figure 1 depicts schematically with the accordingly subscripts. Sh,Eh,IhandRhamass the number of people at each stage of the model description, and Sv,Ev,Ivsigni es proportion of vectors (0Sv; Ev; Iv1 andSv+Ev+Iv= 1). Demographic changes in the number of humans are not considered because the timescale of infection is much faster than the timescale of birth and deaths for the latter to signi cantly alter the development of the disease in the Brazilian context (supplementary material B provides tools to test this assertion). The total vector population is maintained constant during the analysis, although variations of the proportion of vectors on the particular compartments are in- troduced via birth and death rates. The vector in question is regarded as a hypothetical mosquito apt to being infected or infectious throughout all its life- time | which means the model accounts only for the adult stage of their life cycle | and also unable to recover. The time elapsed while an individual is on the aforementioned exposed group is known as the latent period of an organism and, in this work, is adopted as equivalent to the commonly called incubation period (the time elapsed between being infected and exhibiting symptoms), since data is extremely sparse on the latter for humans [35] (namely, the intrinsic incubation period). Both terms are used interchangeably hereafter, and the concepts do not di er on the mosquitoes 3 SvShEhIhRh Iv Evγ δNv δδ αvαh δHuman Population Vector Population(βvIh)/Nh(βhIv)/NvFigure 1: Schematic representation of the SEIR{SEI model for the Zika virus outbreak de- scription. case (the extrinsic incubation period) [36]. In addition, all the members of the susceptible group are treated as equally capable of being infected and the recovered ones as completely immunized. 2.2. Model equations The evolution of individuals through the SEIR-SEI groups is governed by the following (nonlinear) autonomous system of ordinary di erential equations dSh dt= hShIv Nv;dSv dt=Nv vSvIh NhSv; dEh dt= hShIv Nv hEh;dEv dt= vSvIh Nh( v+)Ev; dIh dt= hEh Ih;dIv dt= vEvIv; dRh dt= Ih;dC dt= hEh;(1) whereNrepresents the total population and 1 = the disease's incubation period (each with the corresponding subscript of hfor human's and vfor vector's), 1 = means the vector lifespan, 1 = is the human infection period | which is de ned in this work as the interval of time that a human is infectious | and identi es the transmission rate, speci cally his the mosquito-to-human rate and vthe human-to-mosquito rate. The transmission terms hShIv=Nvand vSvIh=Nhare composed by a number of susceptible individuals ( Sh; Sv, respectively), a transmission rate ( h; v) and the probability of the contact being made with an infectious mem- ber of the other population ( Iv=Nv; Ih=Nh). Both transmissions terms come 4 from the assumption that the rate of contacts is constant, which characterizes a frequency-dependent transmission [37]. The compartmentalization hypothesis requires setting the variables at the initial time of the analysis, t0, such that their sum equals the total population in each case, e.g, Sv(t0) +Ev(t0) +Iv(t0) =Nv= 1. By so doing, for all t>t 0, Sh(t),Eh(t),Ih(t) andRh(t) will always add up to Nh; andNv(t) will be the equilibrium solution of the initial value problem dNv dt=(1Nv(t)); N v(t0) = 1; (2) namely,Nv(t) = 1. This consideration allows the simpli cation of treating the total vector population as a parameter Nv, instead of a variable, since it stays constant throughout the analysis. If one wishes to treat the vector population as a variable, a recruitment parameter would need to be added in place of Nv as well as another di erential equation to account for the changes in Nv(t). The dC=dtequation allows evaluation of the cumulative number C(t) of infectious people until the time t, that is, the amount of humans so far that contracted the disease and have passed through or are in the infectious group at the given time. Additionally, a set of M= 52 points to represent the number of new infec- tious cases of Zika fever at each week is de ned as follows: Nw=CwCw1;N1=C1; w = 2;;52; (3) whereCwis the cumulative number of infectious humans in the w-th EW. Figure 2 organizes the data of cumulative number of infectious and new cases per week provided by the Brazilian Ministry of Health [38] (supplemen- tary material A) for 2016, where the evolution of the infection can be seen. The expected behavior for the model is that it generates a seemingly year long 1020304050 time (weeks) 0 100 200 300number of people103 (a) Cumulative number 0 6 13 19 25number of people103 1020304050 time (weeks)data (b) New cases Figure 2: Outbreak data. outbreak for 2016 that matches C(t) andNwto the data of Figure 2. The end 5 of this outbreak is marked by the exposed and infectious portions of vectors reaching very low values, e ectively killing o the contamination process to the point thatSh(t) stabilizes. 2.3. Nominal parameters The preliminary values for the parameters of the set of Eqs. (1) come from the related literature concerning the Zika infection, vector-borne epidemic mod- els, the Aedes aegypti mosquito (which is the main vector for Zika, Dengue and Yellow fever in Brazil) and publications provided by health organizations and government agencies. The Brazilian Institute of Geography and Statistics re- ports that Brazil had approximately Nh= 206106people by July, 2016 [39], andNvis set 1 to entail an entire vector population. The adopted extrinsic in- cubation period is 1 = v= 9:1 days [25]. This value agrees with other statistical con dence intervals (CI) that are presented for the parameter in another works (95% CI: 7:3{9:3 days [40]) and is close to the numbers suggested by experimen- tal studies for the time necessary for the virus to reach the mosquito's saliva after an infectious blood meal (5 [41, 42] and 7 days [43]). A systematic review and pooled analysis of the literature and case studies available in [44] estimates that the median intrinsic incubation period is 5 :9 days (95% CI: 4 :4|7:6). This values is selected for 1 = hin this work and is compatible with the range of 3{12 days recommended by multiple sources [15, 45, 46], also formerly used in pre- vious studies [47]. The aforementioned literature analysis in [44] also concludes that 9:9 (95% CI: 6 :9-{21:4) days is the mean time until an infected has no detectable virus in blood. Considering the assumption that the infectiousness in Zika infection ends 1 :5{2 days before the virus becomes undetectable [25, 40], the chosen value for the human infectious period is 1 = = 9:92 = 7:9 days. As for the vector lifespan 1 =, \the adult stage of the mosquito is considered to last an average of eleven days in the urban environment " according to [48]. This is the assumed value for the parameter in this work, which is also consistent with the usual life expectancy for the mosquito in Rio de Janeiro, Brazil [49], and comes close to the average of 2{3 weeks considered in biological studies about the species [50] and by the Centers for Disease Control and Prevention [51]. Lastly, the time between a mosquito being infected and it infecting a human, 1= h, and the time between a human infection and a mosquito taking an infec- tious blood meal, 1 = v, is estimated in [40] as an average of 11 :3 days (95% CI: 8:0{16:3) and 8:6 days (95% CI: 6 :2-11:6), respectively. 3. Calibration of the epidemic model 3.1. Forward Problem The epidemic model of section 2, supplemented by an appropriate set of initial conditions, is a continuous-time dynamical system _x(t) =f x(t);p ;x(t0) =x0 (4) 6 where x(t) = Sh(t);Eh(t);Ih(t);Rh(t);Sv(t);Ev(t);Iv(t);C(t) 2R8is the vec- tor of states at time t,x0= Si h;Ei h;Ii h;Ri h;Si v;Ei v;Ii v;Ci 2R8is a prescribed initial condition vector referring to the initial time t0of the analysis, the vec- torp= (Nh; h; h; ;N v; v; v;)2R8lumps the model parameters and f:UR8R8!R8is a nonlinear map which gives the evolution law of this dynamics, de ned (for xed t) on the open set U=n (x(t);p)2R8R8 xn(t)>0 andpn>0;forn= 1;;8o :(5) The forward problem consists in providing initial conditions (IC) and a set of parameters, represented by the pair = (x0;p), and compute by means of numerical integration the model response x(t) from which a scalar observable ( ;t) is obtained. In the forward problem, represents all IC and system parameters from Eq.(1), while ( ;t) is the new cases Nwsystem response from Eq.(3). Since the map fhas a polynomial nature in x, it is a continuously di eren- tiable function in x. Thus, the existence and uniqueness theorem for ordinary di erential equations guarantees that the initial value problem of (4) has an unique solution. Besides that, one can also show that this solution is continu- ously dependent on , as well as the forward map [52, 53]. The evaluation of the system response in the forward problem is performed numerically in this work via a Runge-Kutta (4,5) method and the scalar observ- able of interestNwis used to assess the simulation when compared with real data of the 2016 outbreak made available by the Brazilian Ministry of Health [38] (supplementary material A). The referred data consists of probable cases of infected people per EW, registered by sanitary outposts and health institutions throughout the country when the common symptoms of Zika fever were exhib- ited by an individual. In accordance with the hypothesis that one only displays symptoms when inside Ih(t), theC(t) variable models this discrete accumulat- ing data on a continuous sense and Nwprovides the corresponding in ux per EW. TheC(t) time series is also observed as a criteria for good tting, since it signi es the general impact of the epidemic and because a reasonable Nwre- sult does not necessarily implies an acceptable cumulative number of infectious individuals for all tin comparison to the data. The initial time of the analysis was established as the rst EW of 2016. The remaining individuals in both populations are assumed susceptible at rst, meaningSi h=NhEi hIi hRi handSi v=NvEi vIi v. The initial values for the exposed and infectious groups are set equal, Ei h=Ii handEi v=Ii v. Likewise, the number of infected humans at the initial time must be Ii h=Ci, given its de nition. The value of Ciis taken as the number of con rmed Zika cases in Brazil on the rst EW of 2016 [38], 8 ;201 individuals, and the recovered group is assumed equal to the suspected number of infected in 2015, according to the data available [13], Ri h= 29;639 individuals. As for the proportion of infectious vectors in the rst week, to work around the lack of data for this initial condition, repetitive manual estimations were tried until the resulted time series ofNwpresented reasonable values compared to the real data. It became clear 7 that the system response is very sensible to Ii v, as slight variations in its value are required to achieve feasible results. In the process of choosing its value, the matching of theNwcurve's peak to the amplitude of infection is also a priority, since this is the main interest region for evaluation of the outbreak. The nominal values of the parameters exhibited viable Nwcurves around Ii v= 2:2104. Figure 3 presents the SEIR-SEI model response for the nominal set of param- eters from section 2.3, supplied with the above IC, on an epidemiological week temporal domain consisting of one to fty-two weeks (7 to 365 days), compared with the data of the outbreak (red dots). 1020304050 time (weeks) 0 100 200 300number of people103 data model (a) Cumulative number 0 6 13 19 25number of people103 1020304050 time (weeks)data model (b) New cases Figure 3: Outbreak data (red) and model response using the nominal parameters (blue). The general shape of the curves in Figure 3 do provide qualitative infor- mation regarding the evolution of the infection, even though the portrayed de- scriptions are not quantitatively realistic. This inherent pattern agreement and numerical mismatch suggests that the model response may di er from the real data due to the use of unsuitable values for the parameters or incorrect IC as- sumptions. The search for parameter and initial condition values that make the simulations t well to the observed data de nes model calibration (or system identi cation), being the object of interest of the next section. 3.2. Inverse Problem The model calibration problem seeks to nd a set of parameters that, to a certain degree, makes the model response as close as possible to the empirical observations (reference data), once, due to the erros on model conception and reference data acquisition, it is (practically) impossible for the forward map to reproduce the outbreak observations. The mathematical setting for this case considers the parameters vector de ned in the parameter space E=R12, since here comprises all IC and system parameters from Eq.(1) , excepting Nh,Nv,Ri handCi, which are kept xed in their nominal values. For the purpose of comparison between observations and predictions, a discrete set with Mtime-instants is consid- ered, so that scalar observations and predictions are respectively lumped into 8 y= (y1;y2;;yM) and ( ) = ( 1; 2;; M), both de ned in the data spaceF=RM. Note that the forward map :E!Fassociates to each pa- rameters vector an observable vector ( ) where the component represents the number of new cases in each week of the year, i.e., w=Nw. In practice, the parameters vector is restricted to be on the convex set of admissible values C=n 2E lb ubo , in which lbandubare lower and upper vector bounds for , respectively. In formal terms, given an observation vector y2Fand a prediction vector ( )2F, the calibration aims at nding a vector of parameters such that = arg min 2CJ( ); (6) for a mis t function J( ) =jjy ( )jj2=8 < :MX m=1 ym m( ) 29 = ;: (7) This is the inverse problem associated to the epidemic model. In general this type of problem is extremely nonlinear, with none or low regularity, multiple solutions (or even none), being much more complicated to attack in comparison with the forward problem [54, 55]. A schematic representation of the forward and the inverse problem associated to the epidemic model is shown in Figure 4. This inverse problem attempts to estimate a nite number of parameters on a nite dimension space, being de ned in terms of a typical nonlinear mis t function. Therefore, Theorem 4.5.1 of Chavent [56] can be invoked to guarantee a proper sense of well-posedness (existence, uniqueness, unimodality and local stability) for the inverse problem. The Trust-Region-Re ective method (TRR) is employed here to numeri- cally approximate a solution for the inverse problem (6). The main idea of the method is to minimize a simpler function that re ects the behavior of J( ) in a neighborhood (trust-region) around . The simpler function is de ned as de- pendent on the trial step s, characterizing the Trust-Region subproblem, and its computation is optimized by restricting the subproblem to a two-dimensional subspace. The subspace is linear spanned by a multiple of the gradient gand (in the bounded case) a vector obtained in a scaled modi ed Newton step, used for the convergence condition D( )2g( ) = 0, where Dis a diagonal matrix that depends on ,g,lb, and ub[57]. Finally, the trial step is found through the subproblem as s= arg min s1 2sTQs+gTsj jjDsjj2 ; (8) where  is a scalar associated with the trust region size; Qis a matrix involving D, a Jacobian matrix (also dependent on ,g,lb, and ub) and an approxi- mation of the Hessian matrix [57]. The quadratic approximation in Eq.(8) has well-behaved solutions [58] and if J( +s)<J( ) then is updated to +s 9 Forward Problem Inverse ProblemMathematical Model ˙x(t) =f(x(t),p) x(t0) =x0IC and Parameters α= (x0,p) α∗ Fitting IC and ParametersModel Observable φ(α,t) y1,y2,...,yM Field ObservationsFigure 4: Schematic representation of the forward and inverse problems associated to the epidemic model. and the process iterates, otherwise  is decreased. In addition, a re ection step also occurs if a given step intersects a bound: the re ected step is equivalent to the original step except in the intersecting dimension, where it assumes the opposite value after re ection. The TRR algorithm also requires an initial guess for each parameter, iden- ti ed in the next section as \TRR input". The stopping criteria are the norm of the step and the change in the value of the objective function, with a tolerance of 107. Supplementary material B provides details on the software used for implementation. 3.3. Numerical experiments for calibration Figure 5 presents the best result for the Nwsystem response tting problem using the nominal parameters and IC from sections 2.3 and 3.1 as initial guesses for the TRR algorithm. The upper and lower bounds used for the parameters were set compatible with the literature suggested intervals and are presented in Table 1, along with the parameters and IC values resulted from the calibration (\TRR output"). The ubforwas assumed lower than the lbfor vto maintain consistency with the model interpretation. The minima for Si handSi vwere set to 0:9Nhand 0:99, respectively, to establish a high number of susceptible individuals as is expected for the beginning of an outbreak. Also, the lower and 10 upper bound (0 :999) forSi vwere motivated by noticing how variations in Ei v andIi vof order 103already bring signi cant changes in the system response. The lack of available data for the exposed and infectious groups at the onset of the epidemic was circumvented by appointing its minimum and maximum possible values as lbandub, i.e.,Ei vandIi vwere restricted between zero and one, whileEi handIi hwere bounded by zero and Nh. Additionally, to ensure the model hypotheses of compartmentalization and constant population, two additional tting points were de ned,  h=Si h+Ei h+ Ii hand  v=Si v+Ei v+Ii v, which were set to match NhRi handNv= 1 on Eq.(6), respectively. However, the algorithm is only capable of approximating hand  vto their intended values. So to account for these minor di erences, the resulting values of ( NhRi h)hand 1vwere added to the TRR output ofSi handSi v. These corrections did not impact the calibration, since the scale of the di erences would always be, correspondingly, 104and 102(at most), which are below the sensibility of Si handSi v; thus, they were only exacted to keep the hypotheses rigorously sustained, otherwise the sum of the compartments in each population would quickly tend to NhandNvin an asymptotic fashion. 1020304050 time (weeks) 0 100 200 300number of people103 data model (a) Cumulative number 0 6 13 19 25number of people103 1020304050 time (weeks)data model (b) New cases Figure 5: Outbreak data (red) and calibrated model response (blue) from Table 1. Clearly, Figure 5 is a reasonable description of the outbreak: the general shape of the infection evolution is attained, all parameters are within realistic possibilities, the C(t) curve overshoots the data by merely 6 :00%, the peak value ofNwdi ers from the empirical data maximum by 7 :87% and is only one week o . However, taking into consideration the order of magnitude of the rst data point ( Ci= 8;201) and the scale of the infection (215 ;319 probable cases until the 52th EW [12]), the TRR output for Ii h(253;360) is probably too high, even though there is no reference value to compare with the number of infectious individuals at the beginning of 2016, making it dicult to ascertain on a deterministic manner what is a feasible value for Ii h. Figure 6 allows examination of the system behavior when the Ii hvalue is aroundCi, by depiction of another result to the inverse problem when the upper bound of the initial number of infectious individuals is set to 10 ;000. The 11 TRR input lb ub TRR output Reference h1=5:9 1 =12 1=3 1=12 [44, 45, 46, 15, 47, 40] v1=9:1 1 =10 1=5 1=10 [41, 43, 40] 1=7:9 1 =8:8 1=3 1=8:8 [47, 25, 40] 1=11 1 =21 1=11 1=16:86 [51, 50, 47, 49] h1=11:3 1=16:3 1=8 1=16:3 [40] v1=8:6 1 =11:6 1=6:2 1=11:6 [40] Si h205;953;959 0:9NhNh 205;700;000 |||{ Ei h8;201 0 Nh 15;089 |||{ Ii h8;201 0 Nh 253;360 |||{ Si v0:99956 0 :99 0:999 1 |||{ Ei v2:21040 1 0 |||{ Ii v2:21040 1 0 |||{ Table 1: TRR setup for the Figure 5 calibrated response. The values referring to parameters are in days1, human IC are expressed in individuals, and vector quantities in proportion. same restriction was made over Ei hmerely to simplify the analysis. Table 2 displays the resulting parameters and IC. The model response for a 10 ;000Ii hrestriction also presents acceptable pre- dictions of the general shape and numbers of the outbreak, even though it is less accurate than Figure 5 on a tting criteria for Nw. For comparison, the Nw peak and data maximum di erence increased to 10 :57% and two weeks, while the overshoot on the C(t) time series actually reduced to 5 :74%. 1020304050 time (weeks) 0 100 200 300number of people103 data model (a) Cumulative number 0 6 13 19 25number of people103 1020304050 time (weeks)data model (b) New cases Figure 6: Outbreak data (red) and calibrated model response (blue) from Table 2. To compare the two systems de ned by the parameters and IC from each table, Figure 7 portrays their Ihresponse. The system from Table 1 has an almost monotonically decreasing Ihcurve, except for a slight local maximum around the 6th EW, while the system from Table 2 presents a signi cant increase 12 TRR input lb ub TRR output Reference h1=5:9 1 =12 1=3 1=12 [44, 45, 46, 15, 47, 40] v1=9:1 1 =10 1=5 1=10 [41, 43, 40] 1=7:9 1 =8:8 1=3 1=3 [47, 25, 40] 1=11 1 =21 1=11 1=21 [51, 50, 47, 49] h1=11:3 1=16:3 1=8 1=10:40 [40] v1=8:6 1 =11:6 1=6:2 1=7:77 [40] Si h205;953;959 0:9NhNh 205;953;534 |||{ Ei h8;201 0 10 ;000 6;827 |||{ Ii h8;201 0 10 ;000 10;000 |||{ Si v0:99956 0 :99 0:999 0:999586 |||{ Ei v2:21040 1 4 :14104|||{ Ii v2:21040 1 0 |||{ Table 2: TRR setup for the Figure 6 calibrated response. The values referring to parameters are in days1, human IC are expressed in individuals, and vector quantities in proportion. in the number of infectious individuals by the same time, which correspond to the weeks right before the peak infection. As stated, the lack of empirical data for the current number of infectious at each EW makes it impossible to determine what is a quantitatively reasonable prediction for Ihvalues. But this work assumes that a more possible scenario involves a Ihtime series that also follows the general shape of an epidemic curve around the weeks of maximum infection, especially considering that Figure 7(a) implicates the notion that most people were infected strictly before 2016, which does not seems the case suggested by the Brazilian outbreak data of probable cases of infected per EW [12] when compared to numbers available for 2015 [13]. Thus, with this qualitative criteria in mind, the system from Table 2 is selected for a further analysis over its behavior dependency to the Ii hinitial condition. 1020304050 time (weeks) 0 100 200 300number of people103 (a)Ih(t) response from Table 1 1020304050 time (weeks) 0 3 7 10number of people103 (b)Ih(t) response from Table 2 Figure 7: Comparison between Ih(t) responses from Table 1 (left) and Table 2 (rigth). 13 Figure 8 displays the C(t),NwandIh(t) responses per EW for various Ii h on the system with parameters from Table 2. To simplify the analysis, Ei his considered equal to Ii hin each case. The remaining IC are the same from the Table. The pattern suggests that a Ii hincrease on the system with this given set of parameters continuously escalates the C(t) andNwcurves, eventually overshooting the data by far, and reduces the variations of Ih(t) curve around its local maximum. Figure 8 allows one to make better predictions about the outbreak by analyzing the multiple possible scenarios of epidemic evolution over di erent values for the IC missing empirical information. Clearly, the system response in all cases is qualitative reliable in simulating the outbreak ( C(t) and Nwshape) and can quantitatively approximate the real data values for some values ofIi h. 1020304050 time (weeks) 0 167 333 500number of people103 Data 1000 500010000 20000 3000040000 50000 60000 (a) Cumulative number 1020304050 time (weeks) 0 10 20 30number of people103 Data 1000 5000 10000 20000 30000 40000 50000 60000 (b) New cases 1020304050 time (weeks) 0 20 40 60number of people103 1000 5000 10000 20000 30000 40000 50000 60000 (c)Ihresponse 1020304050 time (weeks) 0 5 10 15number of people103 1000 5000 10000 20000 30000 40000 50000 60000 (d)Ihresponse magni ed Figure 8: Multiple C(t),NwandIh(t) responses per EW using the parameters from Table 2. The values on the legend correspond to the used Ii handEi hto generate the curves. The other IC are also from Table 2. The bottom right graph is the Ihcomparison magni ed around the local maximum region. The red circles are the previously used real data. 14 4. Concluding remarks A SEIR-SEI epidemic model to describe the dynamics of the 2016 Zika virus outbreak in Brazl is developed and calibrated in this work. Nominal quantities for the parameters are selected from the related literature concerning the Zika infection, the Aedes aegypti genus of mosquitoes, vector-born epidemic mod- els and information provided by health organizations. The calibration process is done through the solution of an inverse problem with the aid of a Trust- Region-Re ective method, used to pick the best parameter values that would t the model response for the number of new infectious cases per week into the disease's empirical data. Results within realistic values for the parameters are presented, stating reasonable descriptions with the curve shape similar to the outbreak evolution and proximity between the estimated peak value and data for maximum number of infected during 2016. Further analysis of the results about the lack of data for an initial condition is performed, exhibiting a range of values over which the system response keeps its quantitative reliability to a certain degree. This work is part of a long project of modeling and prediction of epidemics related to the Zika virus in the Brazilian context [23, 24]. In upcoming studies the authors intend to take into account the uncertainties underlying the model parameters via Bayesian updating and employ an Active Subspace approach [59, 60, 61] to explore relevant scenarios in parametric studies. Acknowledgments The authors are indebted to the Brazilian agencies CNPq (National Council for Scienti c and Technological Development), CAPES (Coordination for the Improvement of Higher Education Personnel) and FAPERJ (Research Support Foundation of the State of Rio de Janeiro) for the nancial support given to this research. They are also grateful to Jo~ ao Peterson and Vin cius Lopes, both engineering students at UERJ, for the collaboration in the initial stages of this work. The anonymous reviewers made a series of comments and suggestions that greatly enriched this paper, for which the authors are very grateful.