This work deals with the calibration of SEIR epidemic model for Zika virus outbreak, in Rio de Janeiro, occurred in 2016, based on sim- ilar numerical analysis applied for disease's evolution in French Polyne- sia islands. Demographical increasing or decreasing is not considered. Vector is assumed to be infected for it's all lifespan. Model parameters tting follows an empirical process, by comparison of obtained results with that in Polynesia case study, featuring as satisfactory the values that best reproduce model evolution in Rio de Janeiro scenario. Keywords : nonlinear dynamics, mathematical biology, SEIR epidemic model, model calibration, Zika outbreak in Rio de Janeiro. 1 Introduction Zika virus was rst isolated in primates from Zika forest at Uganda about 1947 [1]. Evidences of the virus in humans were reported in Nigeria at 1968 [2]. It's rst observation in south America occurred in Easter island, about 2014 [3]. First cases in Rio de Janeiro, Brazil, were reported in 2015 [4], quickly evolving to an outbreak in 2016. This scenario caused concern in medical community, health authorities and local population, specially due to a probable relation among Zika virus and other diseases, such as microcephaly [5, 6] and Guillain-Barr e syndrome [7], and the realization of the 2016 Summer Olympics at Rio de Janeiro, which, in reason of the large ux of people, could spread Zika virus to several parts of the world. 1eber.paiva@uerj.br 2michel.tosin@uerj.br 3joao.peterson@uerj.br 4vinicius.lopes@uerj.br 5americo@ime.uerj.br 2 In this (near) epidemic scenario, developing strategies of control and pre- vention of the disease is a critical issue. In this sense, a mathematical model able to predict the number of infected people, during the virus outbreak, is an useful tool, that can be used to identify vulnerable aspects of the disease control strategies. This paper deals with the study of Zika virus outbreak at Rio de Janeiro in 2016. For this purpose, a mathematical model that was used to describe the outbreak in French Polynesia [8] is adapted to Rio de Janeiro scenario. The nominal value of model parameters are estimated based on real data about the outbreak. Numerical simulation results are compared to the out- break ocial numbers, in order to calibrate and validate the adapted model. The rest of this paper is organized as follows. In section 2, the math- ematical model is described, as well as its parameters. In section 3, the process of estimation of model parameters is discussed and numerical sim- ulation results are reported, being followed by the comparison of numerical predictions and experimental data in section 4. Finally, in section 5, the main contributions of this work are emphasized and paths for future works are suggested. 2 Mathematical model The transmission of Zika virus occurs, mainly, by Aedes genus of mosquitoes, but recently researches indicates that Culex genus is able to act as vector too [9]. Medical literature also relates sexual contact and blood transfusions as potential ways of transmission [10]. But, for purposes of modelling, in this work only the transmission by Aedes mosquitoes is considered. To describe the Zika virus outbreak in Rio de Janeiro, this work employs an epidemic model of SEIR (susceptible, exposed, infected and recovered) type [11], adapted from the one used in French Polynesia [8]. In this model, the total population is distributed in four groups (susceptible, exposed, in- fected and recovered) according to the health condition of the individual at timet. The number of susceptible individuals at time tis denoted by S(t), and comprehends healthy individuals which does not have any immunologi- cal resistance against Zika virus. The label of exposed individuals, denoted byE(t), applies to the cases which are in an incubation period (latent pe- riod), i.e. a time gap while, despite of pathogen exposure, individuals are not infectious (asymptomatic period). The number of infected individuals, dubbedI(t), means a ected but not latent anymore. This group is, how- ever, capable of transmitting the disease. The recovered or immune group, 3 SvSh Eh Ih Rh Iv Evγ δ δδ(βvIh)1N αvβhIv αh δHuman Population Vector Population Figure 1: Schematic representation of SEIR epidemic model used to describe Zika virus outbreak. represented by R(t), contains the individuals that were exposed and are recovered. The evolution of the diseased vector population obeys almost the same rules, in exception that there is no recovered group and it's lifetime is consid- ered, since it is small compared to human lifetime. The individuals stream trough the model groups is depicted in Figure 2, both for a ected and vector populations. The stream (evolution) of individuals trough the groups follows the rules de ned by the following set of di erential equations dSh dt= hShIv; (1) dEh dt= hShIv hEh; (2) dIh dt= hEh Ih; (3) dRh dt= Ih; (4) dC dt= hEh; (5) 4 dSv dt= vSvIh NSv; (6) dEv dt= vSvIh N(+ v)Ev; (7) dIv dt= vEvIv; (8) where represents the transmission rates and , the recovery rate, de ned as the inverse of infection period; means disease's incubation ratio, calcu- lated as the inverse of incubation period and , the vector lifespan, and N represents human total population. The hsubscripts refers to the human parameters, while v, the mosquitoes (vector). 3 Numerical experiments 3.1 Model parameters The following values are used, at rst, for the model parameters: N= 16:5106humans, a recent estimation for Rio de Janeiro population [12]; extrinsic incubation ratio v= 1=14 days1[13]; intrinsic incubation ra- tio h= 1=7:5 days1[14]; human infectious ratio h= 1=7 days1[15]; mosquito lifespan ratio v= 1=25 days1[16]; while the vector-to-human transmission rate h= 400 days1, and human-to-vector transmission rate v= 40 days1were arbitrarily chosen through model experimentation, starting from the Tahiti values proposed for the French Polynesia out- break [8]. Regarding the initial conditions, a fully susceptible population is con- sidered, meaning Sh(0) =Nfor humans and Sv(0) = 1 for the proportion of mosquitoes. The hypothesis made for the French-Polynesian model, that initial number of infected (humans or vector) was equal to their respective incubating group is adopted: Ih(0) =Eh(0) andIv(0) =Ev(0) [8]. As such, the value of Iv(0) = 0:014 from the Tahiti outbreak is also used [8]. The initial number of infected humans considered is Ih(0) = 412, such number is a sum of sparse reported cases of Zika infection through some Rio de Janeiro cities in 2015 [17]. Of course, their sum is not a accurately representation of the quantity Ihat such time (beginning of 2016), but a arbitrary choice to condensate the 2015 information due to the lack of more precise data. 5 Lastly,R(0) is arbitrarily chosen to be zero, signifying no recovered indi- vidual at the beginning of the outbreak, and C(0) =Ih(0) because of the system dynamics implications. 3.2 Model predictions In order to verify the representativeness of the proposed model, it is nec- essary to compare its predictions with eld data. In this sense, real data about Zika outbreak, extracted from the epidemiological reports of Brazil's Ministry of Health [18] is used as reference. The report catalog data about the zika infection in the form of suspected accumulated cases since the thir- teenth epidemiological week (03/27/2016 until 04/02/2016). The precise number of con rmed cases is actually unknown by the health authorities. Thus, this work took this value as 35%, as a rst estimate. Unfortunately, only six points of signi cant data are available through these reports so far (depicted on Table 1), and Zika cases along with the associated neurological and congenital syndrome related were not given no- ti able conditions by the Brazil Ministry of Health until February 17th of 2016 [19], thus leaving 2015 data unreliable. week number of cases 13 25930 16 32312 18 38196 20 43516 21 46027 32 60176 Table 1: Number of suspected cases of infected humans in 2016 [18]. The comparison between the number of new infected humans per week, predicted by the model, and the real data is shown in Figure 2. Clearly, the reader can see that model response (blue curve) is not in good agreement with the real data (red circles). Being 2015 the year of the ocial rst con rmed case of Zika in Brazil and the probable epidemic peek that happened in mid-July of the same year [17], the lack of information about the spread of the disease during that year proposed a major diculty on judging the consistency of the model. The choosing of parameters by ways of physical and biological context was not enough to establish the system on a basis of reasonable comparison with empirical data, and the lack of such data also hindered proper selection of 6 010 20 30 40 50010002000300040005000 time (weeks)number of new cases/weekmodel validation Figure 2: Comparison between the model prediction for the number of new infected humans per week (blue curve), and real data (red circles) obtained from [18]. parameters as h, vand initial conditions as IV(0). Therefore, a calibration process is required to analyze the necessary parameters to t the system. 4 Model calibration To calibrate the model, a non-linear curve tting procedure, based on a trust-region-re ective algorithm [20], is employed. The algorithm is per- formed rst varying a single parameter, while the others remained xed. This procedure allows to identify the parameters that most a ect the model response. The parameters hand vshow up to be unreliable, since they required very high initial guesses and usually would deform the quality of the tting curve. Besides that, the model is very sensible to variations on hand v, being more a ected by the values of h. The other parameters do not bring remarkable conclusions. Next, a two varying parameters attempt is conducted through the method. Most combinations do not bring satisfactory results, to the extent that some pairs of parameters could not even be computed in the adopted tolerance levels, e.g ( h; v), probably because of inecient initial guesses or the low number of data points available. Furthermore, it sounds that some parame- ters control the quality of the tting curve. Such control parameters are h and vand, unsurprisingly, the pair ( h; v) proved to be the best one for tting purposes. Hence, the pair ( h; v) is chosen to be focus of the process of model calibration. The parameter Nis held xed in all attempts. A change in the initial 7 conditions is also performed, considering in particular IhandIv. The ini- tial number of infected humans proved unresponsive, threatening signi cant changes on the curve only when varying his order of magnitude, what is beyond of our assumptions. But Ivis noteworthy. The value of Iv= 0:008 is chosen when testing the ( h; v) pair due to compatible general results. Some considerations on this process of investigation to nd a good tting curve regard the model curve peak position (local maxima). Brazil's 2016 Zika peak infection happened on February [19], but precise data for the state of Rio de Janeiro can only be estimated thus far, as discussed. This implicates that the six data points of Table 1 give very few information on the general shape the curve should have around the peek, because of their distance from such event. Additionally, the lack of reliable data about the infection in the beginning of year 2016 and end of 2015 allows a multitude of di erent shape-wise curves to t the midyear data provided by the Ministry of Health. These problems manifested themselves while trying the parameter testing method, as a great number of low peak or wrong peak-positioned curves would t the precise empirical numbers for the region (Table 1) but, of course, contradict the other information about the general distribution of the infection along the year. The adopted solution for this lack of information is to perform (educated) guesses for the the peak region values, around the month of February. An imprecise number of 3500 new cases per week on February for Rio de Janeiro that circulated around brazilian media [21] is used on the iterative algorithm in search of a hand vthat would better represent a peak like curve around such month. Figure 3.a summarizes the result of only using the data from Table 1, where the results of the estimation algorithm are h= 0:0017 days1and v= 4:0599 days1. Figure 3.b represents the curve when a single unocial data point (marked by a asterisk) of 3500 new cases per week is assumed on the epidemiological week 7. For this case the tting algorithm returned in h= 0:0055 days1and v= 0:6479 days1. Finally, Figure 3.c depicts the best tted curved with respect to the empirical data at disposal, considering three extra unocial data points: in addition to the 3500 new cases per week on the seventh week, a 2000 new cases/week on the second and 1300 weekly new cases on the third are included. The resulting estimated parameters are h= 0:0045 days1and v= 0:8296 days1. Figure 3.a indicates a general shape of the curve by only considering the data of Table 1: a early peak on the year 2016 and a decay like ap- pearance around the data points. The introduction of an arbitrary point in Figure 3.b moved the peak closer to the desired week and de ned better the shape. However, without comparison data around the maximum region, the 8 algorithm could not direct the parameters to values that would make 3500 be about the highest new infected per week around the seventh week. Using three arbitrary points as disposed on Figure 3.c, intentionally chosen to do the shape of the peak, the curve better distributed itself through the data and the magnitude stayed closer to the expected maximum value. 010 20 30 40 50010002000300040005000 time(weeks)number of new cases/weekmodel validation (a) 010 20 30 40 50010002000300040005000 time(weeks)number of new cases/weekmodel validation (b) 010 20 30 40 50010002000300040005000 time(weeks)number of new cases/weekmodel validation (c) Figure 3: Camparision between the model prediction (blue curve) and data from epidemiological bulletin (red circles). In the top only bulletin data is presented, while in the bottom 'ghost data' (black stars) are also used. Although graphical results above are satisfactory, the estimated param- eters are not good representations of reality. A h= 0:0039 translates to a vector only contacting a human being in a period of 256 days, which sur- pass extensively the lifespan of the mosquito. This may indicate this SEIR model is unrealistic to represent Rio de Janeiro outbreak, once it ts the observation data only when unreal parameters are used. 9 5 Final remarks This work adapted a SEIR epidemic model, used to predict French Polyne- sia Zika spread, to describe the 2016 Zika virus outbreak in Rio de Janeiro, Brazil. Nominal quantities for the parameters were select through available experimental data in conjunction with arbitrary adopted values. A model calibration procedure, that uses a nonlinear curve tting method, was em- ployed to pick the best parameter values that would t the model response into the empirical data regarding the evolution of the infection. Due to the low number of real data about the Zika outbreak, the curve tting proce- dure was not e ective at rst. The e ects of the lack of information have been solved by incorporating the reference data ansatz for the outbreak peak values. In this second step, the calibration procedure could get a good t to describe the outbreak. Nevertheless, the adjust parameters obtained this way are unrealistic, which may suggest that this SEIR model is not a realistic representation of the Rio de Janeiro outbreak scenario. In future work the authors intend to construct a probabilistic model of uncertainties, to better describe the model parameters variability, and employ a Bayesian procedure for model calibration, in order to to verify if the model calibration can only be done via unrealistic data. Acknowledgments The authors are indebted to the Brazilian agencies CNPq, CAPES, and FAPERJ for the nancial support given to this research.