MODELING EPIDEMIC OUTBREAKS IN GEOGRAPHICAL REGIONS: SEASONAL INFLUENZA IN PUERTO RICO

. We develop a model for the spatial spread of epidemic outbreak in a geographical region. The goal is to understand how spatial heterogeneity inﬂuences the transmission dynamics of the susceptible and infected populations. The model consists of a system of partial diﬀerential equations, which indirectly describes the disease transmission caused by the disease pathogen. The model is compared to data for the seasonal inﬂuenza epidemics in Puerto Rico for 2015-2016.


1.
Introduction. Mathematical modeling of the spatial spread of epidemic diseases in geographical regions presents many challenges. Many authors have developed continuum partial differential equation models to describe spatial variation in epidemics, and typically, diffusion terms are adopted to describe the spatial movement of susceptible and infected individuals in the epidemic populations (see e.g. [1,4,5,7,8,11,12,14,20,21,22,23,24,25,26]). The main difficulty of this approach is that human movement in space does not correspond to random motion assumed by continuum diffusion formulations.
A recent treatment of this problem was given in [16], in which infected individuals diffused in space, but susceptible individuals did not. The underlying assumption was that the diffusion of infected individuals indirectly accounted for the spatial propagation of the microbial agent causing the infection. This assumption is reasonable for outbreak diseases such as seasonal influenza, in which infectiousness of infected individuals lasts a short time compared to the duration of the epidemic, and results primarily in the infection of nearby susceptible individuals. The model in [16] was applied to the seasonal influenza epidemics in Puerto Rico in 2016-2017, and gave reasonable agreement with reported case data for these epidemics.
In this work we develop another approach to account for spatial behavior in outbreak epidemics. We assume that neither susceptible nor infected individuals diffuse in the given spatial region, but infectious individuals infect nearby susceptible individuals, as described by a spatially dependent transmission term. In this way the underlying microbial agent, which is not modeled directly, spreads spatially in the geographical region. Our model also allows for an incubation period before infected individuals become infectious, which is modeled by a time delay in the spatially dependent transmission term.
Yearly outbreaks of seasonal influenza lead to 3-5 million infected cases and approximately 375,000 deaths per year world-wide [28]. Despite the tremendous research effort on seasonal influenza, it is still not well-understood how the seasonal pattern of influenza emerges [19]. We refer readers to the survey paper [13] for proposed theories on the causes of the seasonality of influenza. We propose here a model of seasonal influenza that suggests that the duration of an influenza epidemic is largely determined by the depletion of the susceptible population to a level that no long sustains the transmission of the infection [16].
Our model will be applied to simulating the spatial spread of the 2015-2016 seasonal influenza in Puerto Rico. We propose one possible explanation for the seasonality of influenza in Puerto Rico, which lacks seasonal distinctions throughout the year. A major influenza-virus subtype is introduced to some location in Puerto Rico each year in the winter (e.g. through bird migration or human importation to Puerto Rico from other places), which then spreads through the island and fades away as it exhausts the susceptible population.
Our main objective in this paper is to demonstrate the possibility of applying spatial epidemic models to the study of outbreak diseases in geographical regions using spatial data, rather than using only total case data independently of spatial locations. We note that for most influenza epidemics, the reported case data is only a small fraction of the total case data [15].
In Section 2 we present the model, which consists of a system of partial differential equations for the susceptible, infectious, and removed infected individuals in the epidemic population. In Section 3, we provide proofs for theoretical results of the model. In Section 4, we apply the model to the seasonal influenza epidemic in Puerto Rico in 2015-2016. A major challenge for this application is to identify the geographical boundaries of Puerto Rico and the spatial discretization of the epidemic populations. In Section 5, we compare the numerical simulation of the model to reported case data. In Section 6, we discuss conclusions and extensions of the model to other epidemics in spatial settings.

Model.
Let Ω ⊂ R 2 be a bounded domain with smooth boundary ∂Ω. Let s(t, x), i(t, x) and r(t, x) be the density of susceptible, infected and recovered individuals at position x ∈ Ω and time t, respectively. We propose the following compartmental model to simulate the spatial spread of influenza in a bounded region Ω: where β > 0 is the disease transmission rate, τ > 0 is the time needed for a patient to become infectious after exposure to the virus, γ > 0 is the disease recovery rate, and h(u) = u p 1 + κu q with p ≥ 1, κ ≥ 0 and q ≥ 1, and b(t, x) is the solution of the elliptic equation where ε > 0 is a positive constant, and ∂b/∂ν is the outward normal derivative on ∂Ω.
Equation (2) expresses that fact the pathogen is spreading around the infectious individuals. The function b(t, x) measures the infectiousness of the infected individuals to nearby susceptible individuals. If i(t, x) = δ(x − x 0 ) is the Dirac measure centered at x 0 ∈ Ω, then ∫ Ω b dx = 1 and b(t, x) describes the capability of an infected individual at position x 0 to infect the individual at x. As → 0, b(t, x) → i(t, x), which means the infectious individual at position (t, x) will only infect an individual at the same position. Remark 1. Equation (2) can be derived from the following parabolic system as η goes to 0: The process of letting η → 0 corresponds to the assumption that the dynamics of the virus is faster than the dynamics of humans. Model (1) coupled with (3) is a variant of the model proposed in [16], which assumed that the susceptible population was not moving, and the transmission of virus was viewed indirectly as the diffusion of the infected population.
We impose the following initial conditions where s 0 , i 0 and r 0 are nonnegative continuous functions. The density of exposed individuals is given by the following formula Then e(t, x) satisfies Adding (1) and (5), we have and we deduce from this formula that where n is the total population and does not change in time. The analysis of the model is provided in section 3. In particular, we prove i(t, x) → 0 as t → ∞, as in [16], which means that the disease will eventually disappear.
3. Analysis of the model. The first result is about the existence of solutions of the model.
Theorem 3.1. Assume that β, γ ∈ C(Ω) are strictly positive, s 0 , i 0 , r 0 are nonnegative continuous functions, and κ, , p, q are positive constants. Then the model (1)-(4) has a unique global nonnegative mild solution (s(t, x), i(t, x), r(t, x)), which is classical for t ≥ τ .  (6), T = ∞ and therefore the solution is global. Moreover by [17] or [18,Theorem 1], the solution is classical for The second result is about the positivity of the solutions of the model.
We next study the global dynamics of the model.  Let (s(t, x), i(t, x), r(t, x)) be the solution of (1)-(4). Then Proof. Adding up the equations for s, e and i, we have Integrating the equation on (0, t), we get for all x ∈Ω. Since 0 ≤ s(t, x), e(t, x), i(t, x) ≤ n(x), ∂ t i(t, x) is bounded by the equation of i. Therefore by (9), i(t, x) → 0 for all x ∈Ω. Integrating (8) over Ω and (0, t) leads to Since ∂ t Ω i(t, x)dx is bounded by the equation of i, we have i(t, ·) → 0 in L 1 (Ω).
Integrating both sides of the first equation of (2), Combining (11)-(12), we have It then follows from (10) that which implies Therefore, r ∞ is bounded and measurable, so it is integrable. By the Lebesgue Theorem, r(t, ·) → r ∞ in L 1 (Ω).

4.1.
Using GPS geographical data for model input.

4.1.1.
Constructing the boundaries of municipalities. For partial differential equations with complex geometries, such as the Puerto Rico geometric region Ω, the finite element method (FEM) is particularly suitable. The island of Puerto Rico consists of 76 municipalities, with a total population of approximatedly 3.5 million peoples, in a geographical region of approximately 170 km by 60 km. The first major difficulty is to construct the boundaries of these 76 municipalities. The GPS coordinates of every municipality were collected from the website [29]. We used the MATLAB PDE tool box, which requires the same points for shared common boundaries of any two municipalities. The MATLAB function used to construct this FEM mesh is pdetool. To increase the computation speed, we reduced the number of points defining boundaries. In order to solve this problem we developed a series of MATLAB procedures to construct the geometric coordinates for the boundaries of municipalities, so that common boundaries are defined by the same sequence of points. In Figure 1 we plot these points.

4.1.2.
Mesh generation and finite element method with MATLAB. We first construct the mesh of the domain Ω with small triangles, which preserve the municipality structure. In order to define the required triangulation, we use the Matlab pdetool. The mesh obtained is plotted in Figure 2. It is important to observe that this mesh preserves the municipalities geometry. We use this mesh to define characteristic functions for each municipality. This mesh will be used to construct the densities of the populations of the municipalities. We briefly describe our numerical scheme for the model. At each nodepoint p in Ω h , we solve the equation (1) satisfied by the pathogen density b(t, x). We use the FEM to obtain the approximate solution B n p at each time t n , at each nodepoint p in Ω h . We then update the approximate values (S n p , I n p R n p ) of (S(t n , p), I(t n , p), R(t n , p)) at time t n+1 = t n + ∆t by the following scheme: where d τ ∆t corresponds to the delay τ . Then B n+1 p is obtain from I n p again by using MATLAB to solve a finite element discretized version of Note that the advantage of (15) is the preservation of the stability of the scheme and the positivity of the solutions without any condition on the time step ∆t.   Figure 3. We will simulate the total influenza cases of the 2015-2016 epidemic season. It is very important to note that the total infected cases are much larger than the total reported cases, as the fraction of unreported cases of influenza is very high.
In Figure 4 we plot the density of infected individuals for week 52 in 2015, which will be used as the initial data i(0, x) (week 1) for the simulation. We also plot the

4.2.2.
Reconstruction of population densities per municipalities. The four major municipalities of Puerto Rico with largest population are the eastern San Juan (population 2,350,000), the southern Ponce (population 262, 000), the western Arecibo (population 193,000), and the western Mayaguez (population 89,000). A continuum density function for n(x) is constructed based on the populations and geographical centers of the 76 municipalities. The density of population in the municipalities is based on the data collected from US Census Bureau.
In Figure 5 we plot the municipalities population densities, by taking into account the area of each municipality. We assume that everybody is susceptible at the beginning of the epidemic season. Therefore, the initial value s(0, x) can be approximated by the density of the population. In fact, it is well-known that a significant fraction of the population is immune to a given influenza strain, and we refer to [15] for consideration of this issue.

5.
Simulation of the epidemic outbreak for the 2015-2016 epidemic season. In this section, we simulate the 2015-2016 seasonal influenza epidemic in Puerto Rico. In Table 1 we list the parameters used for the simulation of the model. Since the incubation period is assumed to be 2 days [3], which is small compared to 1 week, we use a constant initial distribution for the simulations.
In Figure 6 and Figure 7 we plot the density of the infected population for the weeks 1, 5, 7, and 11 of the 2015-2016 seasonal inlfuenza epidemic. In Figure 8 we observe that the influenza epidemic spreads around San Juan in the interior region of Puerto Rico. We also observe almost no spreading around Ponce in the south of the island. One limitation of this simulation is that we do not observe during this period any significant spreading of the epidemic to the west side of the island, while some spreading to this region is observed in the data.
In Figure 9 the epidemic seams to persist longer on the west side of the island, while the number of cases decreases in the central part of the island. We observe that for this simulation, the spreading of the epidemic is not given by a single    front, but by multiple wave fronts. In this sense, the simulations are similar to Departamento de Salud de Puerto Rico data. 6. Conclusion and discussion. The spatial progression of an epidemic outbreak in a spatially heterogeneous population presents substantial technical difficulties  for mathematical modeling. Diffusion formulations in partial differential equation models have been used by many authors to describe the movement of people in social geographical settings. Such formulations, however, are far from reality, since in most contexts, people do not diffuse in their spatial environments in the same way that particles diffuse in their spatial environments. Alternative formulations of the movement of people, such as network models, patch models, agent based models, and other approaches also have incomplete connection to the way people move.
In this paper, we propose a spatial epidemic outbreak model, which uses diffusion in a different way. In our approach, diffusion does not model the movement of people, but rather models indirectly the transmission of a disease causing microbial agent from infected individual to uninfected individual. The infectiousness of an infected individual is stipulated by a transmission term b(t, x) that describes the infection capacity of disease transmission, depending on the local density of infected individuals and susceptible individuals at time t and position x. We thus focus on the spread of infection, rather than the movement of people. We analyze the dynamics of the model, and prove that the density of the infected population i(t, x) converges to zero, which means that the epidemic will eventually disappear. This result is in agreement with the result for the diffusive SIR model in [16].
Our model is applied to the simulation of the 2015-2016 seasonal influenza epidemic in Puerto Rico. We use geographical and population data for the simulation, and we describe the epidemic dynamics in each municipality of the island of Puerto Rico. Our simulation compares with the case data reported by Departamento de Salud de Puerto Rico.
Based on our studies, we propose one possible explanation for the seasonality of influenza in this example: a subtype of influenza is introduced to some locations of Puerto Rico annually in the late fall or winter; the infection then quickly spreads through the island and exhausts the susceptible population to a level that no longer sustains the transmission of the disease by the late winter or early spring.
We notice that there is a late second peak in total reported cases graphed in Figure 3, and we believe that spatial heterogeneous distribution of the population is one possible explanation for it. Another possible explanation for this small second peak is that an outbreak of type B strain at week 21 in 2016 may account for it (see Figure 10).
Our studies show that it is possible to model spatial epidemics with comparison to geographical, population, and epidemic data. Our simulations are meant to connect theoretical spatial epidemic models with real epidemic outbreak data. At this stage, there are still many difficulties to overcome in order to make spatial epidemic models more realistic. For example, there are many unreported cases for influenza, so it is necessary to understand better the role of unreported cases [15]. There is also an incomplete understanding of the fraction of the population that is susceptible, since many people have immunity to a current influenza strain through previous infection, vaccination, or other reasons. In future work we will address these and other issues concerning the spatial progression of outbreak epidemics. Figure 10. The total number of reported cases of influenza strain subtypes in 2015-2016. An outbreak of type B strain peaks at week 21 in 2016, which may account for the small second peak in total reported cases graphed in Figure 6.