An outbreak vector-host epidemic model with spatial structure: the 2015–2016 Zika outbreak in Rio De Janeiro

Background A deterministic model is developed for the spatial spread of an epidemic disease in a geographical setting. The disease is borne by vectors to susceptible hosts through criss-cross dynamics. The model is focused on an outbreak that arises from a small number of infected hosts imported into a subregion of the geographical setting. The goal is to understand how spatial heterogeneity of the vector and host populations influences the dynamics of the outbreak, in both the geographical spread and the final size of the epidemic. Methods Partial differential equations are formulated to describe the spatial interaction of the hosts and vectors. The partial differential equations have reaction-diffusion terms to describe the criss-cross interactions of hosts and vectors. The partial differential equations of the model are analyzed and proven to be well-posed. A local basic reproduction number for the epidemic is analyzed. Results The epidemic outcomes of the model are correlated to the spatially dependent parameters and initial conditions of the model. The partial differential equations of the model are adapted to seasonality of the vector population, and applied to the 2015–2016 Zika seasonal outbreak in Rio de Janeiro Municipality in Brazil. Conclusions The results for the model simulations of the 2015–2016 Zika seasonal outbreak in Rio de Janeiro Municipality indicate that the spatial distribution and final size of the epidemic at the end of the season are strongly dependent on the location and magnitude of local outbreaks at the beginning of the season. The application of the model to the Rio de Janeiro Municipality Zika 2015–2016 outbreak is limited by incompleteness of the epidemic data and by uncertainties in the parametric assumptions of the model.

Zika virus appears to be primarily spread through the human population through bites from Aedes mosquitos. The virus incubates in a human host over an asymptomatic period lasting from three to twelve days and once fully developed, the virus disease persists for about a week. It is characterized by low grade fever, rash, joint pain, and conjunctivitis (red eyes). Typically it is mild and seldom requires hospitalization. However the virus has two severe complications which make it a menace to public health. The virus has been linked to an increased risk of Guillian-Barre syndrome which is a severe autoimmune disorder [3]. Perhaps even more serious is its linkage to microcephaly birth defects in newborn babies [4].
Zika epidemics are both year-round and seasonal, dependent upon the year-round prevalence or seasonality of the resident mosquito populations. A recent study [5] describes in detail the potential spread of Zika epidemics into African and Asian-Pacific regions by the importation of infected people. The generation of Zika epidemics by the importation of infected people into year-round or seasonal environments is a major public health concern. Recent mathematical models have been developed to understand these concerns [6][7][8][9][10][11][12][13]. We develop a model that describes both year-round and seasonal hostvector epidemic population dynamics in a geographical region. The disease is borne by vectors to susceptible hosts through criss-cross dynamics in a region of spatially distributed vectors and hosts. The epidemic outbreak begins with the arrival of a small number of viremic hosts in one or more locales in which the disease is not yet present. Our goal is to aid understanding of how the introduction of a small number of infected hosts, in a specific location in a geographic region, will result in a dissipated or a sustained epidemic. The focus of the study is examine the influence of spatial effects on these possible outcomes.
We formulate a criss-cross reaction-diffusion partial differential equations model to describe the spatial evolution of an epidemic. Criss-cross reaction-diffusion models for the circulation of disease between vectors and hosts have been used to describe the spatial spread of malaria [14], the spatial spread of Dengue outbreaks [15,16], and the spatial spread of other diseases by many authors [17][18][19][20][21][22][23][24]. We apply our model to the 2015-2016 Zika seasonal outbreak in the urban area of Rio de Janeiro Municipality in Brazil. We numerically simulate the model to analyze varied scenarios of Zika seasonal epidemics in Rio de Janeiro, dependent upon the input of local spatial outbreaks at the beginning of the season and the time-limitation of seasonality.

Methods
The geographical region is denoted by ⊂ R 2 . The background population of uninfected and susceptible hosts in has geographic density H u (x, y), which is assumed unchanging in time in the demographic and epidemic context of the outbreak. Thus, the model is viewed as applicable to an early phase of the epidemic, during which the epidemic does not alter the local geographic and demographic population structure of hosts. The model consists of the following compartments: The density of infected hosts H i (t, x, y) at time t at (x, y) ∈ , with initial condition H i0 (x, y). The density of uninfected vectors V u (t, x, y) at time t at (x, y) ∈ , with initial condition V u0 (x, y).
The density of infected vectors V i (t, x, y) at time t at (x, y) ∈ , with initial condition V i0 (x, y).

Equations of the model
The equations of the model in the case that transmission from vectors to hosts is year- In addition, the following boundary and initial conditions are satisfied: The spatially dependent parameters of the model are as follows: λ(x, y) is the loss rate of the infected host population (due to recovery or other removal). β(x, y) is the breeding rate of the vector population. μ(x, y) is the loss rate of the vector population due to environmental crowding. σ 1 (x, y) is the transmission rate of uninfected hosts and σ 2 (x, y) is the transmission rate of uninfected vectors. The transmission terms for both hosts and vectors are assumed to be of density-dependent form, rather than frequency-dependent form [25]. A comparison of the two forms for spatially dependent models is given in [26]. Since we assume, during the early phase of the epidemic, that the populations of infected hosts (infected vectors) are relatively small fractions of the populations of uninfected hosts (uninfected vectors), the two forms are essentially the same. δ 1 (x, y) and δ 2 (x, y) are the diffusion rates of the infected hosts and infected vector populations, respectively. In the Appendix we prove the well-posedness of the model.

The local basic reproduction number
Define the local basic reproduction number of the model (1), (2), (3) as follows: R 0 (x, y) is interpreted as the average number of new cases generated by a single case at a given location (x, y) in . An analysis of local reproduction numbers for spatially dependent models is given in [27] and in [28]. Our motivation for this definition is the basic reproduction number R 0 of the spatially independent model (Appendix). Simulations of the spatially dependent model show the following behavior: (1) If R 0 (x, y) < 1 everywhere in , then the populations of both infected hosts and infected vectors extinguish, and the populations converge to the disease free equilibrium. (2) If R 0 (x, y) > 1 in some subregion from an initial local outbreak to an endemic equilibrium in , even if the average value of R 0 (x, y) in all of is < 1.

Equations of the model when the vector population is seasonal
If the vector population is seasonal, then equations of the model must be modified to account for seasonality. We assume that the vector population breeding term β(t, x, y) is dependent on time. We assume that in addition to the vector loss parameter μ(x, y) corresponding to carrying capacity, there is a time-independent vector loss term μ 1 (x, y), corresponding to the average vector life-span 1/μ 1 (x, y). The modified equations are  Brazilian Health Ministry [31] reported that Rio de Janeiro State (population approximately 16,000,000) registered a count of 60,176 cumulative cases from January 1, 2016 to August 13, 2016 (incidence of approximately 364 cases per 100,000 inhabitants). In [32] the weekly case data for Rio de Janeiro Municipality is given from November 1, 2015 through April 10, 2016, during which time the reporting of cases became mandatory. The cumulative number of reported cases in the Municipality during this period was 25,400 [32] (incidence of approximately 423 cases per 100,000 inhabitants).

Parameterization of the Rio de Janeiro model
We simulate the model (1), (4), (5) for Rio de Janeiro Municipality with some parameters assumed. The available epidemic data used for comparison to our simulations for the Rio de Janeiro Municipality 2015-2016 Zika outbreak is very limited. Further, the number of unreported cases, necessarily unknown, is a limitation of the applicability of the model for this application. A more precise fitting of parameters μ, σ , and β requires much higher data accuracy specific to the Zika epidemic in the Municipality. Our purpose is to provide a qualitative description of a typical vector-borne epidemic spatial outbreak, and our simulation of this particular outbreak, with its limitations on parameterization, serves this purpose.
Explanations for our assumptions on specific parameter values are as follows: The time units for our simulations are weeks. The spatial units are kilometers and = (−25, 25) × (−12, 12). The boundary conditions for are a reasonable simplification of the costal boundaries and the less populated northern boundary of Rio de Janeiro Municipality. The average length of the infectious period of infected people is approximately 1 to 2 weeks and we set λ(x, y) = 1.0 [33,34]. The average lifespan of female Aedes aegypti mosquitoes is approximately two weeks in an urban environment [35,36], and we set μ 1 (x, y) = 0.5. The total uninfected host population is approximately approximately 6,000,000, with geographical density function H u (x, y) = 50.0+10 2 (1.0+sin(0.02π x) cos(0.03π y)) ( Fig. 2a), which corresponds approximately to the population density distribution in Fig. 1.
We set the density dependent mosquito loss function μ(x, y) = 0.0015(1.0 + 100 gauss(20.0, 30.0, x) × gauss(0.0, 30.0, y)) ( Fig. 2b), which corresponds to higher levels of mosquito control in the eastern region of the Municipality, where the population density is highest. Here gauss(m, sd, x) is the probability density function in x Fig. 2 a The population of susceptible people H u (x, y) in Rio de Janeiro Municipality, which agrees approximately with the geographical population density in Fig. 1. b The spatially dependent mosquito loss function μ(x, y), which is higher in locations of higher population density due to mosquito control measures of the normal distribution function with mean m and standard deviation sd. Set the transmission parameters σ 1 (x, y) = 0.00000049, σ 2 (x, y) = 0.78 (we assume that individual mosquitoes bite multiple people, people receive multiple bites, and the probability of infection of mosquitoes is much higher than the probability of infection of people).
The diffusion terms for the infected people, uninfected mosquitoes, and infected mosquitoes in the model are understood as idealizations of the indirect spatial spread of the Zika virus infection agent. The spatial spread of the virus is dependent on the direct spread of infected people and uninfected/infected mosquitoes. The spatial movement of people in an urban setting is extremely complex, and a major challenge for epidemic modeling. We set the infected people diffusion parameter δ 1 = 0.2, which provides a simplified way of describing the movement of infected people, in the context of the epidemic, with respect to the spatial spread of the virus. We set the mosquito diffusion parameter δ 2 = 0.2, which is consistent with an estimated adult mosquito dispersal of 30 − 50 m per day [36].
For simplicity, we assume that the mosquito life-span is independent of spatial location, and also independent of time in the season, although for some Aedes species, in some environments, the life-span is correlated to temperature [35]. We take the time dependent mosquito breeding function as β(t, x, y) = 300.0 emg(t,μ,σ ,λ), where emg is the shifted exponentially modified gaussian Here Erfc is the complementary error function. The parameters areμ = −2.0,σ = 5.0, λ = 0.2. The graph of the seasonal mosquito breeding function β(t) is given is Fig. 3.
The assumptions on the parameters of the mosquito population yield a very rapidly rising population at the beginning of the season, which quickly stabilizes to maximal capacity of approximately 14 million, and then declines gradually to very low levels from midway through the season to the end of the season. The total mosquito population, both uninfected and infected, remains mostly uniformly spatially distributed throughout the Municipality throughout the season. The infected mosquito population spatial distribution is very similar to the spatial distribution of infected people. During much of mosquito season, the ratio of total mosquitoes to total people is approximately 2 to 1, which agrees with the ratio in [15].
We set the initial outbreaks in variable locations in the Municipality. For the initial spatial distribution of infected people we set H i (0, x, y) = H i0 gauss(x 0 , 1.0, x) × gauss(y 0 , 1.0, y), centered at (x 0 , y 0 ). The initial number of infected people at the location (determined by H i0 ) is viewed as small and above a threshold level capable of outbreak. It includes imported cases (first order) and possibly some cases generated by first order cases (higher order).

Simulations of the model for Rio de Janeiro
We provide four simulations of the model with initial outbreaks in different locations in the Municipality. The simulation agrees qualitatively with the weekly reported case data for Rio de Janeiro Municipality in [32] (Fig. 4). The spatial distribution of infected people expands from a very small number of initial cases in a small eastern subregion of the Municipality, and disperses throughout the eastern region of the Municipality (Fig. 5). In Fig. 6 we graph the total number of infected people and the total cumulative number of infected people throughout the season. The simulation agrees qualitatively with the weekly reported case data for Rio de Janeiro Municipality given in [32], with approximately 25,500 reported cases between week 44, 2015 and week 15, 2016.
Example 2. We repeat the simulation with the only change from Example 1 the location of the initial outbreak. We take the initial outbreak location as the center of the Municipality with x 0 = 0 and y 0 = 0, R 0 (0, 0) ≈ 1.26. The total number of infected people at time 0 is 20 (H i0 = 20). The infected population again expands from the initial location and disperses throughout the eastern region of the Municipality, but at approximately one-tenth of the number of infected cases as in Example 1 (Figs. 7 and 8). The reason is that R 0 (x, y) is much lower in this initial location than the initial location in Example 1, and the rise of the epidemic is much slower than in Example 1. In Fig. 8 (bottom) we repeat the example with the center of the outbreak location at x 0 = −10, y 0 = 0, R 0 (−10, 0) ≈ 0.52. The infected cases decrease rapidly to 0, because R 0 (x, y) is even lower in the region of the outbreak.
Example 3. We again repeat the simulation with the only change from Example 1 the location of the initial outbreak. We take the initial outbreak with two locations in the center of the Municipality with   The total number of infected people at time 0 is 30. The infected population again expands from the initial location and disperses throughout the eastern region of the Municipality, but at approximately one-third of the number of infected cases as in Example 1 (Figs. 9 and 10). The reason is that the R 0 (x, y) is again lower in the initial outbreak locations than the initial outbreak location in Example 1.
Example 4. We again repeat the simulation with the only change from Example 1 the location of the initial outbreak. We take the initial outbreak with three locations in the eastern region of the Municipality with For initial data in the eastern region (where R 0 (x, y) > 1), the number of infected cases increases and converges to an endemic steady state. For initial data in the western region (where R 0 (x, y) < 1), the number of infected cases first decreases, and then increases to the same endemic state. The simulations indicate the importance of spatial heterogeneity in epidemic models, especially for outbreak scenarios. The importation of a small number of infected cases to isolated localities, may at first dissipate in sub-regions with R 0 (x, y) < 1 , but later rise and spread to sub-regions with R 0 (x, y) > 1, and establish endemicity in the greater geographical region.

Discussion and conclusions
The model (1), (2), (3) describes criss-cross vector-host transmission dynamics of an epidemic outbreak in a geographical region , where the vector population is present year-round. The outbreak occurs with a small number of infected hosts in a small subregion of the much larger geographical region . The diffusion terms describe the on-going average spatial spread of the disease microbial agent within infected vectors and infected hosts in the geographical region. The focus of the model is to describe the geographical spread from an initial localized immigration into the region, in terms of the epidemiological properties of the outbreak vector-host transmission dynamics.
We prove that the partial differential equations model (1), (2), (3) is mathematically well-posed, and compare its properties to an analogous ordinary differential equations model in the spatially independent case (Appendix). The outcomes of the model depend on the spatially distributed local reproduction number R 0 (x, y). In the case of year-round vector settings, simulations indicate that the connection of R 0 (x, y) to the outcome of an outbreak is as follows: if R 0 (x, y) < 1 everywhere in , then the epidemic will extinguish; if R 0 (x, y) > 1 in some subregion of , then the epidemic has the possibility to spread from an initial outbreak to an endemic equilibrium in , even if the average value of R 0 (x, y) < 1 throughout all of .
The model Eqs. with the initial data located in the eastern region. Bottom: Model simulation of the spatial distributions of infected cases at time=0, 20, 40, 130, with the initial data located in the western region. Both simulations converge to the same limiting density, but the one with the initial data in the western region first decreases before increasing and converging with the reported case data in the Municipality [32]. We argue that the assumption of an unchanging number for the susceptible population is reasonable for the Zika outbreak in Rio de Janeiro Municipality. The justification for this assumption is based on current demographic data for the Municipality [37]. Between 2010 and 2016 the population increased from approximately 6.000,000 at approximately 0.49% per year. The total number of reported cases during the 2015-2016 outbreak is less than 1% of the susceptible population, which is not significantly depleted during the outbreak.
A limitation of our model is the difficulty of estimating the number of unreported cases, and in some examples of Zika epidemics the ratio of reported cases to unreported cases has been quite high. In one study, the Federated States of Micronesia in 2007, the number of reported cases was 108 and the number of unreported cases (estimated through seroconversion testing) was estimated at 74% of the total population of 7,391 [38]. In another study, the French Polynesia outbreak in 2013-2014, the number of reported cases was estimated at 7-17% of the total number of infections, with 94% of the total population infected [33]. The setting for Rio de Janeiro Municipality is very different, however, and the demographic changes in Rio de Janeiro Municipality in one year could off-set a relatively higher ratio of unreported-to-reported cases, given that the reported cases  R 0 (x, y). The average value of R 0 (x, y) over the whole spatial region is ≈ 0.984. Right: The cumulative total number of infected cases in the whole region as a function of time, which for both initial data in the eastern region (blue) and initial data in the western region (red) converge to ≈ 1, 630 represented approximately 0.4% of the population [31,32]. Additionally, the probability of Zika re-infection is not yet fully known. Whether Zika could become established as an endemic disease in a larger urban population thus remains unclear [33]. Our model simulations are based on the number of reported cases, but we note that if the ratio of unreported to reported cases is significantly higher, then the parameters must be adjusted.
A limitation of our model is that it does not take into account the possibility of sexual transmission of Zika. It is noted in [2], however, that sexual transmission is a small percentage of total transmission, and may not initiate or sustain an outbreak. Another limitation of our model is that we assume the uninfected mosquito population is uniformly geographically distributed at the beginning of the season, since there is no detailed temporal geographic mosquito data available for Rio de Janeiro Municipality. We note that current investigations are developing such data for geographical regions, which could be implemented eventually for spatial models of vector borne epidemics as described by our model. One such investigation is Project Premonition [39], developed by Microsoft to autonomously locate, robotically collect, and computationally analyze mosquito populations for pathogenicity in geographical environmental regions.
The model simulation suggests that the Zika epidemic in Rio de Janeiro Municipality may rise each season from initial outbreak locations, with very small numbers of infected people, and spread through a larger region of the Municipality. Although the epidemic subsides at the end of the season, the final size of the epidemic at the end of the season depends on the initial outbreak locations of infected cases in the region, when geographic heterogeneity and time-limited seasonality are taken into account. The local reproduction number R 0 (x, y) indicates that the most effective interventions decrease the infection rates σ 1 (x, y), σ 2 (x, y), increase the isolation of infected people λ(x, y), increase the mosquito removal rate μ(x, y), and control the importation of infected people, all concentrated in regions of high density population H u (x, y) and in the beginning of the season.
For the Zika epidemic in Rio de Janeiro Municipality the model suggests that the outbreak in the 2015-2016 season will occur again in the 2016-2017 season, and in future seasons. The importation of infected cases into the Municipality at the beginning of the season is inevitable, because of the general influx of people into this major metropolitan center of Brazil. Some of these cases will not generate a further spread of cases, but some will, with consideration of spatially variable factors. The reduction of future, and more extensive, seasonal outbreaks of Zika in the Municipality requires higher level monitoring of the people arriving in the region and higher level mosquito control measures throughout the region, again with consideration of spatially variable factors. and initial conditions Proof. We first observe that a unique classical solution {H i (t), V u (t), V i (t)} exists in C 1 ( ) on a maximal interval of existence [ 0, T max ) [40][41][42]. Standard arguments [42] guarantee that {H i (t), V u (t), V i (t)} remain nonnegative for t ∈[ 0, T max ). Moreover, the classical solution can be globally defined if we can establish uniform a priori bounds. Set M(t, x, y) = V u (t, x, y) + V i (t, x, y) and add Eqs. (2) and (3) to obtain Theorem 1 in [24] guarantees the existence of a unique global classical solution M(t) ∈ Further, in [24] it is proved that there exists M ∈ C 0 + ( ), M = 0, such that lim t→∞ M(t) = M ∈ C 0 + ( ). We note that the disease free equilibrium of (1), (2), (3) is (0, M, 0). From [24] there exists N 1 > 0 such that max t≥0 M(t) C 0 Consequently, the solution exists globally on [ 0, ∞).