Soumaya Naanani

Jean-Philippe Monchoux

Catherine Mabru

Alain Couret
email: alain.couret@cemes.fr

Pure climb of [001] dislocations in TiAl at 850 °C

Keywords:

This paper presents a study by Transmission Electron Miscroscopy of the deformation microstructure in the γ phase of a TiAl alloy crept at 850 °C under 150 MPa. A never observed population of dislocations is evidenced and characterized. It is shown that their Burgers vector is of [001] type and that they are moving by pure climb in the (001) planes. The reasons of the presence of these dislocations are discussed.

Although TiAl alloys are candidate to be used at temperatures higher than 800 °C in aeronautic and automotive engines, the deformation mechanisms activated in this temperature range in industrial alloys with complex microstructures are still not really understood. The few studies reported in literature on the deformation of TiAl alloys at high temperatures all agree to attribute the deformation mainly to twinning and ordinary dislocations [START_REF] Appel | Gamma Titanium Aluminides: Science and Technology[END_REF][START_REF] Kad | [END_REF] and more scarcely to [001] and 〈100] dislocations [3][4][5] and 〈011] and 1/2 b 112] superlattice dislocations [6]. For the case of ordinary dislocations, glide configurations characterized by dislocation elongation along their screw direction and by numerous pinning points are still observed at 900 °C [START_REF] Kad | [END_REF]. In aluminum rich single crystals, Inui et al. have noted the glide of these ordinary dislocations in non-compact planes, as (001) and (110) planes [7]. However, a contribution of climb has been evidenced in several contributions [START_REF] Kad | [END_REF]7,8] and associated to the transition between the brittle to ductile material behavior [START_REF] Appel | Gamma Titanium Aluminides: Science and Technology[END_REF][START_REF] Kad | [END_REF]5,8]. Interestingly, Kad and Fraser have shown that a decrease of the strain rate promotes the activation of climb instead of glide [START_REF] Kad | [END_REF]. They have proposed that at low strain rates, climb is activated in a second step when an appropriated density of dislocations has been previously created by glide. Following the simultaneous observations of glide and climb dislocations in samples crept at 700 °C, the activation in some grains of a high density of gliding dislocations similar to that observed at intermediate temperatures has been ascribed to the accommodation of high local internal stresses [9]. From in situ heating experiments, climbing ordinary dislocations are seen to be preferentially oriented along the 〈110〉 orientations and forming helicoidal structure [8,10,11]. Malaplate et al. [10] has demonstrated that during creep at 700 °C, these ordinary dislocations can move by a mixed climb mechanism which is a mixture of glide and climb and which may explain the high measured values of the activation volume, previously underlined by Appel et al. [8]. The [001] dislocations have been observed to glide in {110] planes in single crystals at 950 °C under favorable orientation conditions [3] and lying in pile-ups at 1000 °C in an aluminum rich alloy [5].

In this context, the present study aims at analyzing the dislocations microstructure and at determining the activated deformation mechanisms in an industrial TiAl alloy exhibiting a complex microstructure and crept at 850 °C. Here, observations and characterizations of a population of dislocations are reported, which, in our knowledge, has never been described in the current literature. The reason of the activation of this new type of dislocations will be discussed in the final section.

The alloy investigated is an IRIS alloy (Ti -48Al -2W -0.08B) obtained by the densification of a pre-alloyed powder by Spark Plasma Sintering [12]. Its microstructure is formed by small lamellar grains surrounded by borders made of single phased γ grains that contain B2 precipitates [13]. The alloy was crept at 850 °C under 150 MPa up to 1.5% of strain. Under this solicitation, the minimum creep rate was 6•10 -8 s -1 and the creep life was 260 h. Thin foils suitable for Transmission Electron Microscopy (TEM) were extracted from this crept sample with their planes cut perpendicular to the tensile axis of the sample. They were examined in a JEOL 2010 TEM using a Tilt-Rotation specimen holder commercialized by the Gatan company, with the aim to obtain observations under various diffracting conditions. This holder allows the following movements of the specimen: (i) a 360°rotation around the direction of the electron beam, which is perpendicular to the thin foil plane at zero tilt and (ii) the usual tilt ±40°around the longitudinal tilt axis of TEM holders, which is situated along the vertical direction on all the micrographs presented in this report. As usual for the γ phase of TiAl alloys, the quadratic L1 0 structure will be treated as an ordered face-centered cubic structure consistently with a c/a ratio close to 1.

Fig. 1 displays the deformation microstructure in a γ grain of a border using the g = (111Þ reflection (g is the vector of the reciprocal lattice normal to the plane used to obtain the image in two beams condition; g is called the diffracting vector in what follows). This microstructure consists of nearly rectilinear dislocations, which appear to be elongated along directions situated at 45°from each other. No pinning points generally observed on gliding ordinary dislocations [14][15][16] are seen on these dislocations. Fig. 1,b shows the stereographic projection of this grain under zero tilt. The perpendicular to the thin foil plane (T), which is also the direction of the loading axis is at 22°of the [001] direction. The deformation microstructure has been studied in ten γ grains of the borders by TEM, and in three of them such kind of rectilinear dislocations has been observed. In Fig. 2, the dislocation loop situated in the red rectangle marked on Fig. 1 is observed under different diffracting conditions. The loop is visible with three reflections of the 〈111〉 type. This indicates that the Burgers vector is neither that of an ordinary dislocation (1/2 b 110] type), nor that of a superdislocation of b101] type.

For the micrographs of Fig. 3, the (001) plane is perpendicular to the electron beam (a) and the g diffracting vectors are all perpendicular to the [001] direction; two of them are of b200] type and the two others of the b220] one. The rectilinear segments are parallel to the projections of the b100] and b110] directions of the crystal. Remarkably, as surrounded red, on each of these micrographs one segment is not visible contrary to the others. Moreover, these invisible segments appear to be parallel to the direction of the corresponding g diffracting vector. This kind of contrast is known to exist for the case of edge dislocations while the g b = 0 extinction condition is fulfilled [START_REF] Edington | Practical Electron Microscopy in Materials Science: 3 Interpretaion of Electron Micrographs[END_REF][START_REF] Williams | Transmission Electron Micoscopy: 3 Imaging[END_REF], as exemplified in the complex metallic alloy ξ -Al -Pd -Mn [START_REF] Feuerbacher | [END_REF]. In this configuration, complete invisibility is only obtained if a second condition extinction g b˄u = 0 is also fulfilled (u is the unit vector along the dislocation line). If it is not, strong residual contrasts are observed, particularly when g b˄u N 0.6 [START_REF] Edington | Practical Electron Microscopy in Materials Science: 3 Interpretaion of Electron Micrographs[END_REF][START_REF] Williams | Transmission Electron Micoscopy: 3 Imaging[END_REF]. This kind of residual contrast results from the bugling by edge dislocations of the diffracting plane, which is perpendicular to the extra half plane and contains the Burgers vector (see Fig. 25.4 of [START_REF] Williams | Transmission Electron Micoscopy: 3 Imaging[END_REF]).

To explain all these TEM observations, the only solution is to have the Burgers vector of this dislocation parallel to the [001] direction (b = [001]), which is perpendicular to all the used diffracting vectors (g b = 0). Considering the situation depicted by the Fig. 3, b (g = [200]), three situations are encountered:

• when u//g//[200], b˄u will be perpendicular to g, the scalar product of g by b˄u is equal to 0 and the dislocation will be invisible. That is the case of the red arrowed segment. • when u//[020] and ⊥ to g = [200], b˄u will be parallel to g, the scalar product of g by b˄u is different of 0 and is in fact equal to ±2. The dislocation segment is visible, due to a strong residual contrast. • when u//[220] or [220], the scalar product of g by b˄u is different of 0 and equal to ± ffiffiffi 2 p . The dislocation is visible.

The three other situations of Fig. 3, c-e can be analyzed in the same way.

As a conclusion, all the TEM experimental observations are explained by considering that the dislocation loop is made of [001] dislocations which are pure edge in character and lying along the b100] and b110] directions.

To determine the plane in which these dislocations are moving, tilt experiments have been performed around various directions. Fig. 4 presents one of them, for which, using the holder rotation, the normal to the (001) plane has been placed perpendicular to the tilt axis. It follows that the trace of this (001) plane is going through this tilt axis and stays in place during the tilt experiment (the trace of the plane is its intersection with the thin foil plane). The dislocation loop is observed under different inclinations, using different image conditions (Fig. 4a to h). For each micrograph, the apparent distance between the two exterior segments is measured in the direction perpendicular to the tilt axis. It is reported as a function of the inclination angle (Fig. 4,i). This curve follows a sinusoidal variation with a maximum situated at +20°. The plane of the loop (marked P in Fig. 4,i) is thus a plane with a trace along the inclination axis and which is normal to the observation plane for a thin foil inclination angle of 70°(90°-20°). In the limit of our experimental uncertainty, it appears that this plane is close to the (001) plane. Therefore, the edge [001] dislocations are moving by pure climb. Their rectilinear aspect suggests that the pure climb mechanism occurs by the nucleation and propagation of pairs of jogs as previously suggested by for the case of ordinary dislocations [8,10]. In work of Jiao et al.

[001], [3] dislocations have been observed edge in character and lying along the [110] direction. The authors have concluded to a glide in the (110) plane, a configuration incompatible with the present situation since in our case the two b110] type planes are nearly perpendicular to the observation plane (Fig. 1b).

In TiAl alloys, the deformation is admitted to be mainly due to ordinary dislocations whatever the temperature and stress ranges [START_REF] Appel | Gamma Titanium Aluminides: Science and Technology[END_REF]. However, it exists in the L1 0 structure only two possible directions for the Burgers vectors of ordinary dislocations which are both situated in the (001) plane. Consequently, this kind of dislocations produces a deformation concentrated in this plane and other deformation systems are necessary to insure a three dimensional deformation, or in other words to fulfill the Von Mises's criterion. Movements of superlattice dislocations are scarcely observed except in single γ-phased alloys obtained with aluminum rich compositions [20][21][22]. Twinning due to the glide of 1/6 b 112] Shockley dislocations in (111) plane is in general activated and brings some deformation out of the (001) plane. The twinning activity requires a nucleation process and the Schockley dislocation propagation. On the one hand, the twin nucleation has been found to be due to the decomposition of Frank dislocations in γ grains [23,[START_REF] Couret | Twinning in Advanced Materials[END_REF] or to other similar dislocation reactions [START_REF] Appel | Gamma Titanium Aluminides: Science and Technology[END_REF], and to occur by the emission of Shockley dislocations at interfaces in lamellar alloy [START_REF] Appel | Gamma Titanium Aluminides: Science and Technology[END_REF]. The second process is easier, which explains why twins are more frequently activated in lamellar alloys. On the other hand, in situ straining experiments have shown that the twin propagation is controlled by the overcoming by the Shockley dislocations of heterogeneous pinning points [16,[START_REF] Messerschmidt | Dislocation Dynamics During Plastic Deformation[END_REF]. Hence, twinning appears to be athermal or at least weakly thermally activated. All these properties explain why twinning is not much activated at 850 °C under 150 MPa in single-phased γ grains. Fig. 4. Tilt experiment. The loop is observed under different inclinations (a to h). As indicated by the double arrows, the distance between two segments is measured in each case. Its variation with the inclination angle follows a sinusoidal variation (i). On the stereographic projection under zero tilt (j), P is the pole of the moving plane of the dislocation.

At this high temperature, diffusion is largely efficient enough to activate climb, as previously observed for ordinary dislocations in particular at low strain rate and low stress [START_REF] Kad | [END_REF]. The small size of the γ grains of the borders (a few microns long and about 1 μm wide - [13]) and the proximity of lamellar interfaces is another factor enhancing the diffusion effects through the presence of numerous interfaces. The climb of dissociated superlattice dislocations with antiphase boundary lying in (001) or (111) type planes must involve local rearrangements of the atoms which complicate the occurring of such a mechanism. Similarly to the present results, at high temperatures in B2 intermetallic alloys (FeAl and NiAl), climb of {001} dislocations has been evidenced instead of movement of 〈111〉 superdislocations [22]. Thus, even if it has never been previously observed in TiAl alloys, a contribution of [001] dislocations seems to be an appropriate way to develop a three dimensional deformation. In addition, the orientation of the studied grain with the (001) axis close to the tensile axis of the creep sample represents a favorable situation. The deformation processes, as the glide of [001] dislocation in {110) planes, or that of b010] in the (001) panes reported in Refs [3,4,5] respectively, might be activated because of the reasons mentioned above under different sample orientations.

In conclusion, in a sample crept at 850 °C under 150 MPa, we have observed a population of pure [001] edge dislocations which are moving in the (001) planes by pure climb. The dislocation segments are rectilinear and elongated along the b100] and b110] directions of this plane, which suggests that climb occurs by the nucleation and lateral propagation of jog pairs, as already evidenced for ordinary dislocations of TiAl alloys. Several factors are favorable for the activation of these dislocations, which were never observed: the high temperature, which promotes diffusion, the low stress that inhibits twinning and the orientation of the grain.

Fig. 1 .

 1 Fig. 1. General view (a) and stereographic projection (b) of the γ grain investigated in the present study. On (b), the poles of the plane are labelled by their direction in the reciprocical lattice, T is the direction of the tensile axis of the crept sample.

Fig. 2 .

 2 Fig. 2. TEM observations of the dislocation loop situated inside the red rectangle of Fig. 1 using three reflections of the 〈111〉 type.

Fig. 3 .

 3 Fig. 3. TEM observations of the dislocation loop situated inside the red rectangle of Fig. 1, using four g diffracting vectors all perpendicular to the [001] direction. The segments out of contrast are surrounded in red.

Acknowledgements

This study has been conducted in the framework of the cooperative projects "IRIS-ANR-09-MAPR-0018-06" supported by the French Agence Nationale de la Recherche (ANR) and "ALTIAERO" project of the "IDEX-ATS-2016-060-CIF-D-DRDV" program supported by the Université Fédérale -Toulouse Midi-Pyrénées, which are acknowledged. The authors want to thank Daniel Caillard for fruitful discussions about TEM contrasts and climbing processes of dislocations.


Eber Dantas

Michel Tosin
email: michel.tosin@uerj.br

Americo Cunha
email: americo@ime.uerj.br

Calibration of a SEIR-SEI epidemic model to describe the Zika virus outbreak in Brazil

Keywords: Zika virus dynamics, nonlinear dynamics, mathematical biology, SEIR epidemic model, model calibration, inverse problem

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The Zika virus is a flavivirus that upon infection in humans causes an illness, known as Zika fever, identified commonly with macular or papular rash, mild fever and arthritis [START_REF] Brasil | Zika virus outbreak in Rio de Janeiro, Brasil: clinical characterization, epidemiological and virological aspects[END_REF][START_REF]Zika Virus[END_REF]. It is mainly a vector-borne disease carried by the genus Aedes of mosquitoes [START_REF]Zika Virus[END_REF][START_REF] Fernandes | Culex quinquefasciatus from Rio de Janeiro is not competent to transmit the local Zika virus[END_REF], while in a lesser amount it is also transmitted via sexual interaction [START_REF] Petersen | Update: Interim Guidance for Preconception Counseling and Prevention of Sexual Transmission of Zika Virus for Persons with Possible Zika Virus Exposure-United States, September 2016[END_REF][START_REF] Coelho | Higher incidence of Zika in adult women than adult men in Rio de Janeiro suggests a significant contribution of sexual transmission from men to women[END_REF], and contamination by blood transfusion is under investigation [START_REF] Motta | Evidence for Transmission of Zika Virus by Platelet Transfusion[END_REF]. The Zika virus was first isolated in primates from the Zika forest in Uganda in 1947 [START_REF] Dick | Isolations and serological specificity[END_REF]. Evidences of the virus in humans were found in Nigeria in 1968 [START_REF] Moore | Arthropod-borne viral infections of man in Nigeria, 1964-1970[END_REF]. An epidemic occurred in 2007 on Micronesia [START_REF] Duffy | Zika virus outbreak on Yap Island, Federated States of Micronesia[END_REF], followed by multiple outbreaks on several Pacific Islands between 2013 and 2014 [START_REF] Cao-Lormeau | French Polynesia, South Pacific[END_REF][START_REF] Tognarelli | A report on the outbreak of Zika virus on Easter Island, South Pacific[END_REF]. The first Zika virus autochthonous case in Brazil was reported around April, 2015 [START_REF]Boletim Epidemiolólico. monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika até a Semana Epidemiológica[END_REF], and nearly 30,000 cases of infection were already notified by January 30, 2016 [START_REF] Faria | Zika virus in the Americas: early epidemiological and genetic findings[END_REF], 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 [START_REF] Hennessey | Zika virus spreads to new areas -Region of the Americas, May 2015-January 2016[END_REF]. The Brazilian Ministry of Health registered 215,319 probable cases of Zika fever (130,701 of which were confirmed) until the 52th epidemiological week (EW) of 2016 [START_REF]Boletim Epidemiolólico. monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika até a Semana Epidemiológica[END_REF]. 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 [START_REF] Petersen | Update: Interim Guidance for Preconception Counseling and Prevention of Sexual Transmission of Zika Virus for Persons with Possible Zika Virus Exposure-United States, September 2016[END_REF][START_REF] Valentine | Zika virus-associated microcephaly and eye lesions in the newborn[END_REF] and Guillain-Barré syndrome [START_REF] Santos | Zika virus and the Guillain-Barré syndrome -Case series from seven countries[END_REF], whose correlation to the Zika virus was considered by the World Health Organization a "scientific consensus" [START_REF]Zika virus, Microcephaly and Guillain-Barré Syndrome -Situation Report[END_REF].

In this epidemic scenario, the development of control and prevention strategies 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 effective and vulnerable aspects on disease control programs [START_REF] Naheed | Numerical study of SARS epidemic model with the inclusion of diffusion in the system[END_REF][START_REF] Huppert | Mathematical modelling and prediction in infectious disease epidemiology[END_REF][START_REF] Lizarralde-Bejarano | Understanding epidemics from mathematical models: Details of the 2010 dengue epidemic in Bello (Antioquia, Colombia)[END_REF]. Furthermore, for an epidemic model to be truly useful it must undergo a judicious process of validation [START_REF] Oberkampf | Verification and Validation in Scientific Computing[END_REF][START_REF] Cunha | Modeling and Quantification of Physical Systems Uncertainties in a Probabilistic Framework[END_REF], which consists in comparing model predictions with real data in order to evaluate if they are realistic. In general, the first 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 first 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 identification and validation of representative models to describe Zika virus outbreaks in Brazil [START_REF] Dantas | A mathematical analysis about Zika virus outbreak in Rio de Janeiro[END_REF][START_REF] Dantas | Zika virus in Brazil: calibration of an epidemic model for 2016 outbreak[END_REF]. For this purpose, a SEIR-SEI mathematical model is adapted to the Brazilian scenario. This specific SEIR-SEI description has been successfully used before for the outbreaks in Micronesia [START_REF] Funk | Comparative analysis of Dengue and Zika outbreaks reveals differences by setting and virus[END_REF] and French Polynesia [START_REF] Kucharski | Transmission dynamics of Zika virus in island populations: a modelling analysis of the 2013-14 French Polynesia outbreak[END_REF]. Some assumptions were also based on similar studies performed over SEIR dynamical systems [START_REF] Wang | Dynamical analysis and perturbation solution of an SEIR epidemic model[END_REF][START_REF] Wanga | Edge-based SEIR dynamics with or without infectious force in latent period on random networks[END_REF][START_REF] Safi | Qualitative analysis of an age-structured SEIR epidemic model with treatment[END_REF][START_REF] Cai | A stochastic SIRS epidemic model with nonlinear incidence rate[END_REF]. 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 comparison to real data of the outbreak as an initial effort 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.

Epidemic model for Zika virus dynamics

Model hypotheses

This work utilizes a variant of the Ross-Macdonald model [START_REF] Smith | MacDonald, and a theory for the dynamics and control of mosquito -Transmitted pathogens[END_REF] for epidemic predictions, separating the populations into a SEIR-SEI framework (susceptible, exposed, infectious, recovered) [START_REF] Brauer | Mathematical Models in Population Biology and Epidemiology[END_REF]33,34]. Each category represents the health condition of an individual inside such group at time t, with respect to the considered 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. S h , E h , I h and R h amass the number of people at each stage of the model description, and

S v , E v , I v signifies proportion of vectors (0 ≤ S v , E v , I v ≤ 1 and S v + E v + I v = 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 significantly 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 introduced via birth and death rates. The vector in question is regarded as a hypothetical mosquito apt to being infected or infectious throughout all its lifetime -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 differ on the mosquitoes case (the extrinsic incubation period) [START_REF] Lessler | Assessing the global threat from Zika virus[END_REF]. In addition, all the members of the susceptible group are treated as equally capable of being infected and the recovered ones as completely immunized.

S v S h E h I h R h I v E v γ δNv δ δ α v α h δ Human Population Vector Population (β vI h )/N h (β h Iv)/Nv

Model equations

The evolution of individuals through the SEIR-SEI groups is governed by the following (nonlinear) autonomous system of ordinary differential equations

dS h dt = -β h S h I v N v , dS v dt = δ N v -β v S v I h N h -δ S v , dE h dt = β h S h I v N v -α h E h , dE v dt = β v S v I h N h -(α v + δ) E v , dI h dt = α h E h -γ I h , dI v dt = α v E v -δ I v , dR h dt = γ I h , dC dt = α h E h , (1)
where N represents the total population and 1/α the disease's incubation period (each with the corresponding subscript of h for human's and v for vector's), 1/δ means the vector lifespan, 1/γ is the human infection period -which is defined in this work as the interval of time that a human is infectious -and β identifies the transmission rate, specifically β h is the mosquito-to-human rate and β v the human-to-mosquito rate. The transmission terms β h S h I v /N v and β v S v I h /N h are composed by a number of susceptible individuals (S h , S v , respectively), a transmission rate (β h , β v ) and the probability of the contact being made with an infectious member of the other population (I v /N v , I h /N h ). Both transmissions terms come 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, t 0 , such that their sum equals the total population in each case, e.g, S v (t 0 ) + E v (t 0 ) + I v (t 0 ) = N v = 1. By so doing, for all t > t 0 , S h (t), E h (t), I h (t) and R h (t) will always add up to N h ; and N v (t) will be the equilibrium solution of the initial value problem

dN v dt = δ (1 -N v (t)) , N v (t 0 ) = 1 , (2)
namely, N v (t) = 1. This consideration allows the simplification of treating the total vector population as a parameter N v , 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 δ N v as well as another differential equation to account for the changes in N v (t).

The dC/dt equation 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 infectious cases of Zika fever at each week is defined as follows:

N w = C w -C w-1 , N 1 = C 1 , w = 2, • • • , 52 , (3)
where C w is 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 [START_REF] Da | Obtenção de número de casos confirmados de zika, por municiṕio e semana epidemioloǵica[END_REF] (supplementary 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 outbreak for 2016 that matches C(t) and N w to the data of Figure 2. The end of this outbreak is marked by the exposed and infectious portions of vectors reaching very low values, effectively killing off the contamination process to the point that S h (t) stabilizes.

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 models, 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 reports that Brazil had approximately N h = 206 × 10 6 people by July, 2016 [39], and N v is set 1 to entail an entire vector population. The adopted extrinsic incubation period is 1/α v = 9.1 days [START_REF] Funk | Comparative analysis of Dengue and Zika outbreaks reveals differences by setting and virus[END_REF]. This value agrees with other statistical confidence intervals (CI) that are presented for the parameter in another works (95% CI: 7.3-9.3 days [START_REF] Ferguson | Countering the Zika epidemic in Latin America[END_REF]) and is close to the numbers suggested by experimental studies for the time necessary for the virus to reach the mosquito's saliva after an infectious blood meal (5 [41,[START_REF] Li | Oral Susceptibility of Singapore Aedes (Stegomyia) aegypti (Linnaeus) to Zika Virus[END_REF] and 7 days [START_REF] Chouin-Carneiro | Differential Susceptibilities of Aedes aegypti and Aedes albopictus from the Americas to Zika Virus[END_REF]). 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/α h in this work and is compatible with the range of 3-12 days recommended by multiple sources [START_REF] Valentine | Zika virus-associated microcephaly and eye lesions in the newborn[END_REF][START_REF] Ioos | Current Zika virus epidemiology and recent epidemics[END_REF][START_REF]Rapid risk assessment: Microcephaly in Brazil potentially linked to the Zika virus epidemic[END_REF], also formerly used in previous studies [START_REF] Villela | Zika in Rio de Janeiro: Assessment of basic reproduction number and comparison with dengue outbreaks[END_REF]. 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 [START_REF] Funk | Comparative analysis of Dengue and Zika outbreaks reveals differences by setting and virus[END_REF][START_REF] Ferguson | Countering the Zika epidemic in Latin America[END_REF], the chosen value for the human infectious period is 1/γ = 9.9 -2 = 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 [START_REF] Otero | A stochastic population dynamics model for Aedes Aegypti: formulation and application to a city with temperate climate[END_REF]. 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 [START_REF] De Oliveira | Daily survival rates and dispersal of Aedes aegypti females in Rio de Janeiro, Brazil[END_REF], and comes close to the average of 2-3 weeks considered in biological studies about the species [START_REF] Nelson | Aedes Aegypti: Biology and Ecology[END_REF] and by the Centers for Disease Control and Prevention [START_REF]Dengue and the Aedes aegypti mosquito[END_REF]. 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 infectious blood meal, 1/β v , is estimated in [START_REF] Ferguson | Countering the Zika epidemic in Latin America[END_REF] as an average of 11.3 days (95% CI: 8.0-16.3) and 8.6 days (95% CI: 6.2-11.6), respectively.

Calibration of the epidemic model

Forward Problem

The epidemic model of section 2, supplemented by an appropriate set of initial conditions, is a continuous-time dynamical system

ẋ(t) = f x(t), p , x(t 0 ) = x 0 (4)
where

x(t) = S h (t), E h (t), I h (t), R h (t), S v (t), E v (t), I v (t), C(t) ∈ R 8 is the vec- tor of states at time t, x 0 = S i h , E i h , I i h , R i h , S i v , E i v , I i v , C i ∈ R 8
is a prescribed initial condition vector referring to the initial time t 0 of the analysis, the vec- 8 is a nonlinear map which gives the evolution law of this dynamics, defined (for fixed t) on the open set

tor p = (N h , β h , α h , γ, N v , β v , α v , δ) ∈ R 8 lumps the model parameters and f : U ⊂ R 8 × R 8 → R
U = (x(t), p) ∈ R 8 × R 8 x n (t) > 0 and p n > 0, for n = 1, • • • , 8 . (5)
The forward problem consists in providing initial conditions (IC) and a set of parameters, represented by the pair α = (x 0 , 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 N w system response from Eq.( 3).

Since the map f has a polynomial nature in x, it is a continuously differentiable function in x. Thus, the existence and uniqueness theorem for ordinary differential equations guarantees that the initial value problem of (4) has an unique solution. Besides that, one can also show that this solution is continuously dependent on α, as well as the forward map φ [52,[START_REF] Perko | Differential Equations and Dynamical Systems[END_REF].

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 observable of interest N w is used to assess the simulation when compared with real data of the 2016 outbreak made available by the Brazilian Ministry of Health [START_REF] Da | Obtenção de número de casos confirmados de zika, por municiṕio e semana epidemioloǵica[END_REF] (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 exhibited by an individual. In accordance with the hypothesis that one only displays symptoms when inside I h (t), the C(t) variable models this discrete accumulating data on a continuous sense and N w provides the corresponding influx per EW. The C(t) time series is also observed as a criteria for good fitting, since it signifies the general impact of the epidemic and because a reasonable N w result does not necessarily implies an acceptable cumulative number of infectious individuals for all t in comparison to the data.

The initial time of the analysis was established as the first EW of 2016. The remaining individuals in both populations are assumed susceptible at first, meaning

S i h = N h -E i h -I i h -R i h and S i v = N v -E i v -I i v .
The initial values for the exposed and infectious groups are set equal,

E i h = I i h and E i v = I i v .
Likewise, the number of infected humans at the initial time must be I i h = C i , given its definition. The value of C i is taken as the number of confirmed Zika cases in Brazil on the first EW of 2016 [START_REF] Da | Obtenção de número de casos confirmados de zika, por municiṕio e semana epidemioloǵica[END_REF], 8,201 individuals, and the recovered group is assumed equal to the suspected number of infected in 2015, according to the data available [START_REF] Faria | Zika virus in the Americas: early epidemiological and genetic findings[END_REF], R i h = 29,639 individuals. As for the proportion of infectious vectors in the first week, to work around the lack of data for this initial condition, repetitive manual estimations were tried until the resulted time series of N w presented reasonable values compared to the real data. It became clear that the system response is very sensible to I i v , as slight variations in its value are required to achieve feasible results. In the process of choosing its value, the matching of the N w curve'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 N w curves around I i v = 2.2 × 10 -4 . Figure 3 presents the SEIR-SEI model response for the nominal set of parameters from section 2.3, supplied with the above IC, on an epidemiological week temporal domain consisting of one to fifty-two weeks (7 to 365 days), compared with the data of the outbreak (red dots). The general shape of the curves in Figure 3 do provide qualitative information regarding the evolution of the infection, even though the portrayed descriptions are not quantitatively realistic. This inherent pattern agreement and numerical mismatch suggests that the model response may differ from the real data due to the use of unsuitable values for the parameters or incorrect IC assumptions. The search for parameter and initial condition values that make the simulations fit well to the observed data defines model calibration (or system identification), being the object of interest of the next section.

Inverse Problem

The model calibration problem seeks to find 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 α defined in the parameter space E = R 12 , since here α comprises all IC and system parameters from Eq.( 1) , excepting N h , N v , R i h and C i , which are kept fixed in their nominal values. For the purpose of comparison between observations and predictions, a discrete set with M time-instants is considered, so that scalar observations and predictions are respectively lumped into y = (y 1 , y 2 , • • • , y M ) and φ (α) = (φ 1 , φ 2 , • • • , φ M ), both defined in the data space F = R M . Note that the forward map φ : E → F associates to each parameters vector α an observable vector φ (α) where the component represents the number of new cases in each week of the year, i.e., φ w = N w . In practice, the parameters vector is restricted to be on the convex set of admissible values C = α ∈ E lb ≤ α ≤ ub , in which lb and ub are lower and upper vector bounds for α, respectively.

In formal terms, given an observation vector y ∈ F and a prediction vector φ (α) ∈ F , the calibration aims at finding a vector of parameters α * such that

α * = arg min α∈C J(α) , (6)
for a misfit function

J(α) = ||y -φ (α) || 2 =    M m=1 y m -φ m (α) 2    . ( 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 [START_REF] Aster | Parameter Estimation and Inverse Problems[END_REF][START_REF] Yaman | A Survey on Inverse Problems for Applied Sciences[END_REF]. 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 finite number of parameters on a finite dimension space, being defined in terms of a typical nonlinear misfit function. Therefore, Theorem 4.5.1 of Chavent [START_REF] Chavent | Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications[END_REF] can be invoked to guarantee a proper sense of well-posedness (existence, uniqueness, unimodality and local stability) for the inverse problem.

The Trust-Region-Reflective method (TRR) is employed here to numerically approximate a solution for the inverse problem [START_REF] Motta | Evidence for Transmission of Zika Virus by Platelet Transfusion[END_REF]. The main idea of the method is to minimize a simpler function that reflects the behavior of J(α) in a neighborhood (trust-region) around α. The simpler function is defined as dependent 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 g and (in the bounded case) a vector obtained in a scaled modified Newton step, used for the convergence condition D(α)

-2 g(α) = 0, where D is a diagonal matrix that depends on α, g, lb, and ub [START_REF] Coleman | An interior, trust region approach for nonlinear minimization subject to bounds[END_REF]. Finally, the trial step is found through the subproblem as

s * = arg min s 1 2 s T Qs + g T s | ||D s|| 2 ≤ ∆ , ( 8
)
where ∆ is a scalar associated with the trust region size; Q is a matrix involving D, a Jacobian matrix (also dependent on α, g, lb, and ub) and an approximation of the Hessian matrix [START_REF] Coleman | An interior, trust region approach for nonlinear minimization subject to bounds[END_REF]. The quadratic approximation in Eq.( 8) has well-behaved solutions [START_REF] Conn | Trust Region Methods[END_REF] and if J(α + s) < J(α) then α is updated to α + s and the process iterates, otherwise ∆ is decreased. In addition, a reflection step also occurs if a given step intersects a bound: the reflected step is equivalent to the original step except in the intersecting dimension, where it assumes the opposite value after reflection.

The TRR algorithm also requires an initial guess for each parameter, identified 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 10 -7 . Supplementary material B provides details on the software used for implementation.

Numerical experiments for calibration

Figure 5 presents the best result for the N w system response fitting 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 ub for δ was assumed lower than the lb for α v to maintain consistency with the model interpretation. The minima for S i h and S i v were set to 0.9 N h and 0.99, respectively, to establish a high number of susceptible individuals as is expected for the beginning of an outbreak. Also, the lower and upper bound (0.999) for S i v were motivated by noticing how variations in E i v and I i v of order 10 -3 already bring significant 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 lb and ub, i.e., E i v and I i v were restricted between zero and one, while E i h and I i h were bounded by zero and N h . Additionally, to ensure the model hypotheses of compartmentalization and constant population, two additional fitting points were defined, 6), respectively. However, the algorithm is only capable of approximating Σ h and Σ v to their intended values. So to account for these minor differences, the resulting values of (N h -R i h )-Σ h and 1-Σ v were added to the TRR output of S i h and S i v . These corrections did not impact the calibration, since the scale of the differences would always be, correspondingly, 10 -4 and 10 -2 (at most), which are below the sensibility of S i h and S i 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 N h and N v in an asymptotic fashion. 1.

Σ h = S i h + E i h + I i h and Σ v = S i v + E i v + I i v , which were set to match N h -R i h and N v = 1 on Eq.(
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 of N w differs from the empirical data maximum by 7.87% and is only one week off. However, taking into consideration the order of magnitude of the first data point (C i = 8,201) and the scale of the infection (215,319 probable cases until the 52th EW [START_REF]Boletim Epidemiolólico. monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika até a Semana Epidemiológica[END_REF]), the TRR output for I i 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 difficult to ascertain on a deterministic manner what is a feasible value for I i h . Figure 6 allows examination of the system behavior when the I i h value is around C i , 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. same restriction was made over E i h merely to simplify the analysis. Table 2 displays the resulting parameters and IC.

---- E i v 2.2 × 10 -4 0 1 0 ---- I i v 2.2 × 10 -4 0 1 0 ----
The model response for a 10,000 I i h restriction also presents acceptable predictions of the general shape and numbers of the outbreak, even though it is less accurate than Figure 5 on a fitting criteria for N w . For comparison, the N w peak and data maximum difference increased to 10.57% and two weeks, while the overshoot on the C(t) time series actually reduced to 5.74%. 2.

To compare the two systems defined by the parameters and IC from each table, Figure 7 portrays their I h response. The system from Table 1 has an almost monotonically decreasing I h curve, except for a slight local maximum around the 6th EW, while the system from Table 2 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 I h values. But this work assumes that a more possible scenario involves a I h time 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 [START_REF]Boletim Epidemiolólico. monitoramento dos casos de dengue, febre de chikungunya e febre pelo vírus Zika até a Semana Epidemiológica[END_REF] when compared to numbers available for 2015 [START_REF] Faria | Zika virus in the Americas: early epidemiological and genetic findings[END_REF]. Thus, with this qualitative criteria in mind, the system from Table 2 is selected for a further analysis over its behavior dependency to the I i h initial condition. 2. To simplify the analysis, E i h is considered equal to I i h in each case. The remaining IC are the same from the Table . The pattern suggests that a I i h increase on the system with this given set of parameters continuously escalates the C(t) and N w curves, eventually overshooting the data by far, and reduces the variations of I h (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 different values for the IC missing empirical information. Clearly, the system response in all cases is qualitative reliable in simulating the outbreak (C(t) and N w shape) and can quantitatively approximate the real data values for some values of I i h . 2. The values on the legend correspond to the used I i h and E i h to generate the curves. The other IC are also from Table 2. The bottom right graph is the I h comparison magnified around the local maximum region. The red circles are the previously used real data.

--- E i v 2.2 × 10 -4 0 1 4.14 × 10 -4 ---- I i v 2.2 × 10 -4 0 1 0 ----

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 models 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-Reflective method, used to pick the best parameter values that would fit 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 [START_REF] Dantas | A mathematical analysis about Zika virus outbreak in Rio de Janeiro[END_REF][START_REF] Dantas | Zika virus in Brazil: calibration of an epidemic model for 2016 outbreak[END_REF]. 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 [START_REF] Constantine | Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies[END_REF][START_REF] Diaz | A modified SEIR model for the spread of Ebola in Western Africa and metrics for resource allocation[END_REF]61] to explore relevant scenarios in parametric studies.

Figure 1 :

 1 Figure 1: Schematic representation of the SEIR-SEI model for the Zika virus outbreak description.

Figure 2 :

 2 Figure 2: Outbreak data.

Figure 3 :

 3 Figure 3: Outbreak data (red) and model response using the nominal parameters (blue).

1

 1 , y 2 , . . . , y M Field Observations

Figure 4 :

 4 Figure 4: Schematic representation of the forward and inverse problems associated to the epidemic model.

Figure 5 :

 5 Figure5: Outbreak data (red) and calibrated model response (blue) from Table1.

Figure 6 :

 6 Figure 6: Outbreak data (red) and calibrated model response (blue) from Table2.

  I h (t) response from Table1

  I h (t) response from

Figure 8

 8 Figure8displays the C(t), N w and I h (t) responses per EW for various I i h on the system with parameters from Table2. To simplify the analysis, E i h is considered equal to I i h in each case. The remaining IC are the same from the Table.The pattern suggests that a I i h increase on the system with this given set of parameters continuously escalates the C(t) and N w curves, eventually overshooting the data by far, and reduces the variations of I h (t) curve around its local maximum. Figure8allows one to make better predictions about the outbreak by analyzing the multiple possible scenarios of epidemic evolution over different values for the IC missing empirical information. Clearly, the system response in all cases is qualitative reliable in simulating the outbreak (C(t) and N w shape) and can quantitatively approximate the real data values for some values of I i h .

Figure 8 :

 8 Figure 8: Multiple C(t), Nw and I h (t) responses per EW using the parameters from Table2. The values on the legend correspond to the used I i h and E i h to generate the curves. The other IC are also from Table2. The bottom right graph is the I h comparison magnified around the local maximum region. The red circles are the previously used real data.

Table 1 :

 1 TRR setup for the Figure5calibrated response. The values referring to parameters are in days -1 , human IC are expressed in individuals, and vector quantities in proportion.

Table 2 :

 2 TRR setup for the Figure6calibrated response. The values referring to parameters are in days -1 , human IC are expressed in individuals, and vector quantities in proportion.

Table 2 Figure 7

 27

: Comparison between I h (t) responses from Table 1 (left) and Table 2 (rigth).

Acknowledgments

The authors are indebted to the Brazilian agencies CNPq (National Council for Scientific 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 financial support given to this research. They are also grateful to João 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.

Supplementary Material

This supplementary material provides working MATLAB scripts for the main document "Calibration of a SEIR-SEI epidemic model to describe the Zika virus outbreak in Brazil." The first section uses "ode45" function to evaluate the forward problem as explained in section 3.1 of the main document. Section 2 solves the inverse problem via a TRR algorithm (main document section 3.2) implemented in the "lsqcurvefit" function of the software. Finally, section 3 exemplifies a variation of the model by introducing vital dynamics for the human population, useful for comparison.

All the codes in this document use the following data set. See supplementary material A. 2 NewCasesData = [ 8 2 0 1 , 8 3 4 4 , 9 5 1 6 , 1 0 4 4 3 , 1 2 9 4 0 , 1 5 1 5 8 , 2 1 2 8 0 , ... 3 2 0 8 5 4 , 2 0 1 3 9 , 1 7 5 2 0 , 1 5 6 4 4 , 1 2 8 6 5 , 1 2 7 1 5 , 1 1 2 9 8 , 1 0 0 3 4 , 8 4 1 0 , ... 4 7 3 9 7 , 6 1 5 4 , 5 0 0 8 , 4 0 8 9 , 2 9 9 7 , 2 9 3 6 , 2 5 1 6 , 1 8 8 8 , 1 5 2 1 , 1 2 8 2 , 1 1 2 6 , ... 5 1 0 1 0 , 9 1 7 , 7 4 1 , 7 0 3 , 6 6 5 , 5 6 4 , 4 5 5 , 4 2 2 , 3 8 0 , 3 8 0 , 2 9 6 , 3 0 5 , 3 1 5 , ...

Forward Problem

The default values for the parameters and initial conditions (IC) are the TRR output from Table 2, section 3.3 of the main document.

22

Inverse Problem

See section 2.3 and 2.4 of the main document for the nominal parameters and initial conditions used as initial guesses here. The default results are those from Table 2, section 3.3 of the main document.

Changing human population

A recruitment term (η N h ) and a mortality rate (µ) are needed to implement a changing human population. The human variables in Eqs. (1) in section 2.2 of the main document become

Also, a new differential equation must be added to the model to track the varying population

and, consequently, its respective initial condition, N h (0), is required for the forward problem. Code 7 provides the right hand side of this variation of the model. An estimation for the values of η and µ (per thousand individuals in a year) can be found in [? ] for 2015 (the most recent officially published mortality and birth rates for Brazil). The solution of the forward problem with this code can be used to observe how the human vital do not affect N w (neither C(t)) in a one year outbreak for this model. Of course, S h (t) will behave different, as expected.