Assessing the spatio-temporal spread of COVID-19 via compartmental models with diffusion in Italy, USA, and Brazil
Malú Grave, Alex Viguerie, Gabriel F. Barros, Alessandro Reali, Alvaro L. G. A. Coutinho
AA SSESSING THE SPATIO - TEMPORAL SPREAD OF
COVID-19
VIACOMPARTMENTAL MODELS WITH DIFFUSION IN I TALY , USA,
AND B RAZIL
A P
REPRINT
Malú Grave
Dept. of Civil EngineeringCOPPE/Federal University of Rio de JaneiroP.O. Box 68506, RJ 21945-970, Rio de Janeiro, Brazil [email protected]
Alex Viguerie
Department of MathematicsGran Sasso Science InstituteViale Francesco Crispi 7, L’Aquila, AQ 67100, Italy [email protected]
Gabriel F. Barros
Dept. of Civil EngineeringCOPPE/Federal University of Rio de JaneiroP.O. Box 68506, RJ 21945-970, Rio de Janeiro, Brazil [email protected]
Alessandro Reali
Dipartimento di Ingegneria Civile ed ArchitetturaUniversità di PaviaVia Ferrata 3, Pavia, PV 27100, Italy [email protected]
Alvaro L. G. A. Coutinho
Dept. of Civil EngineeringCOPPE/Federal University of Rio de JaneiroP.O. Box 68506, RJ 21945-970, Rio de Janeiro, Brazil [email protected]
February 16, 2021 A BSTRACT
The outbreak of COVID-19 in 2020 has led to a surge in interest in the mathematical modelingof infectious diseases. Such models are usually defined as compartmental models, in which thepopulation under study is divided into compartments based on qualitative characteristics, with differentassumptions about the nature and rate of transfer across compartments. Though most commonlyformulated as ordinary differential equation (ODE) models, in which the compartments depend onlyon time, recent works have also focused on partial differential equation (PDE) models, incorporatingthe variation of an epidemic in space. Such research on PDE models within a Susceptible, Infected,Exposed, Recovered, and Deceased (SEIRD) framework has led to promising results in reproducingCOVID-19 contagion dynamics. In this paper, we assess the robustness of this modeling frameworkby considering different geometries over more extended periods than in other similar studies . Wefirst validate our code by reproducing previously shown results for Lombardy, Italy. We then focuson the U.S. state of Georgia and on the Brazilian state of Rio de Janeiro, one of the most impactedareas in the world. Our results show good agreement with real-world epidemiological data in bothtime and space for all regions across major areas and across three different continents, suggestingthat the modeling approach is both valid and robust. K eywords COVID-19 · Compartmental models · Diffusion-reaction · Partial differential equations · Adaptive mesh refinement and coarsening a r X i v : . [ q - b i o . P E ] F e b PREPRINT - F
EBRUARY
16, 2021
The COVID-19 pandemic has caused an extraordinary global disruption, both in terms of human lives and economicdamage. To limit the spread of contagion, governments have taken unprecedented measures. Such actions have included,but not been limited to, quarantines, curfews, lockdowns, and domestic and international travel suspensions. Whileundoubtedly effective and considered necessary by many experts, these measures also carry a high cost in terms ofeconomic impact and mental and physical health, among others. In part, such restrictions are motivated by the lack ofreliable data and models on this disease’s transmission in time and space, forcing cautious (and not always entirelyrational) responses from the authorities and population. More than ever, these events demonstrate the need for toolsto predict the spatio-temporal dynamics of contagion. It is important to highlight that such tools should be not onlyaccurate, but, given the difference in the quality of available data in distinct regions of the world (or even within thesame country), they should also be robust with respect to the data used for parameter identification. In this context,model validation - also against data possibly coming from different countries - should play an important role.The mathematical modeling of epidemics is a well-established field, and has been applied with success for manyimportant cases. Examples of past applications include Dengue [1], Zika [2], HIV [3], SARS [4], Malaria [5], Ebola [6],and many others. However, the urgency of the COVID-19 pandemic has motivated the need for more research in thisarea, with several models for this pandemic outbreak being presented in recent months [7, 8, 9, 10, 11, 12, 13, 14, 15, 16].Such models typically follow a compartmental framework, in which the population under study is divided intocompartments with assumptions about the nature and time rate of transfer from one compartment to another [17, 8].The large majority of the proposed compartmental models are composed of a system of ordinary differential equations(ODEs) in time. Though such ODE models are simple to formulate, analyze, and solve numerically, they do notnaturally account for individuals’ movements from one region to another. To incorporate spatial variations, one maydefine regional compartments corresponding to different geographic areas and introduce coupling terms to accountfor movements across regions. Such approaches have also been applied to COVID-19 in, for example, [18, 10, 19].ODE models may be used in combination with other approaches as well, such as with a neural network [20] or anagent-based computational framework [21, 22, 23]. While all these approaches are effective, they lack the ability tomodel spatial transmission at the continuous level, which may have some advantages.To this end, one may instead resort to a compartment model based on partial differential equations (PDEs). Suchmodels are decidedly less common than ODE ones, in particular, due to their increased difficulty and time involvedin both implementation and numerical solution. However, such an additional cost is counterbalanced by the factthat PDE models naturally (and accurately) incorporate spatial information, allowing for a continuous descriptionof spatio-temporal dynamics, with the potential to account for geographical features, population heterogeneity, andmulti-scale dynamics with relative ease. While quite uncommon, PDE models have been already used with success inthe past to model epidemics [24, 25, 26].For the specific case of COVID-19, PDE models were applied in [7, 8, 27, 14], in which a compartmental SEIRD model (susceptible, exposed, infected, recovered, deceased) was used. A modified version of this model, with the addition ofa quarantined compartment, was applied in [16]. Such modeling framework has shown promising results, being able torecreate observed data and make predictions with reasonable accuracy, but these successful applications have been sofar obtained on a limited scale, with results restricted to single, specific regions. Therefore, despite the clearly highpotential shown by this modeling approach, important questions remain regarding its robustness . In particular, it is stillto be proven the model capability to guarantee reliable results for a wide range of areas, which may exhibit significantlydifferent geographical, behavioral, and demographic characteristics. Furthermore, even if the model may demonstratethe ability to represent different regions accurately, it is unclear to what extent parameters must be specifically tuned toobtain accurate results for a given case. In fact, ideally, a robust model should be applicable to a wide range of settings,requiring relatively little specific parameter tuning to produce acceptable results for a given case.Within this context, the current work aims at evaluating the robustness of the model used in [7, 8, 27] by studyingin detail its application to a series of different cases. In order to accomplish this task, the model is applied to thefollowing three cases: The Italian region of Lombardy, the U.S. state of Georgia, and the Brazilian state of Rio deJaneiro. These are not only the homelands of the authors, but they do represent regions of three distinct continents withvery different features in terms of culture, population, and mobility, as well as government response to the pandemic. Todeeply test its robustness to different region-specific characteristics, the model is applied across each case with limitedparameter tuning and case-specific adjustments, and it is shown to be capable of producing consistently reliable results,in particular for predictions relative to short to mid-time spans, which are relevant for timely political decisions.The remainder of this work is organized as follows: Section 2 introduces the governing equations and the relevantparameters, as well as the adopted terminology and notation. This section is also used to briefly discuss considerationsregarding the numerical solution of the associated discrete problem. In the subsequent sections, the model is then tested2
PREPRINT - F
EBRUARY
16, 2021on the cases of Lombardy, Georgia, and Rio de Janeiro, and for each case extensive discussion and analysis of resultsare provided. Finally, the paper is concluded by a summary of the main findings and suggestions for future researchdirections.
We present the governing equations following the continuum mechanics framework first shown in [7], as opposedto more traditional notations common in mathematical and biological references. We consider a system that may bedecomposed into N distinct populations: u ( x , t ) , u ( x , t ) , ..., u N ( x , t ) . Let Ω ∈ R be a simply connected domain ofinterest and denote its boundary as ∂ Ω = Γ D ∩ Γ N . Let [0 , T ] denote a generic time interval. The transient nonlineardiffusion-reaction system of equations written in compact vector notation then reads as follows: ∂ u ∂t + ( A + B ( u )) u − ∇ · ( ν ∇ u ) − f = 0 in Ω × [0 , T ] (1) u = u D in Γ D × [0 , T ] (2) ( ν ∇ u ) · n = h in Γ N × [0 , T ] , (3)where s ( x , t ) , e ( x , t ) , i ( x , t ) , r ( x , t ) , and d ( x , t ) denote the densities of the susceptible , exposed , infected , recovered ,and deceased populations, respectively. We additionally denote the cumulative number of infected as c ( x , t ) and thesum of the living population as n ( x , t ) ; i.e., n ( x , t ) = s ( x , t ) + e ( x , t ) + i ( x , t ) + r ( x , t ) . We then use the compactvector notation to represent these quantities as: u = [ s, e, i, r, d ] T .The tensors A , B and ν , and the vector f are defined according to the particular dynamics of the system in question.In most applications, ν = ν ( x , t ) , that is, diffusion is time-dependent, heterogeneous, and anisotropic; however,this is not strictly necessary. In addition to the boundary conditions (2), (3), we must also specify an initial condition u ( x ,
0) = u . We define the total population U i ( t ) of a given compartment u i ( x , t ) as: U i ( t ) = (cid:90) Ω u i ( x , t ) d Ω (4)for i = 1 , , · · · , N .To adequately define a model of COVID-19 contagion, we make the following assumptions, as discussed in [7, 27]:• We consider only mortality due to COVID-19;• We do not consider new births;• We assume that a portion of exposed individuals will never develop symptoms (asymptomatic cases), andhence will move directly from the exposed compartment to the recovered compartment;• We assume that both pre-symptomatic and asymptomatic (the exposed compartment), as well as symptomatic(the infected compartment) individuals may spread the disease;• There is a latency period after exposure and before the development of symptoms;• Movement is proportional to the population size; i.e., we expect greater diffusion to occur in heavily populatedregions;• There is no movement among the deceased population.Under the above assumptions, the frequency-dependent system of equations reads: ∂s∂t + β i n si + β e n se − ∇ · ( nν s ∇ s ) = 0 (5) ∂e∂t − β i n si − β e n se + ( α + γ e ) e − ∇ · ( nν e ∇ e ) = 0 (6)3 PREPRINT - F
EBRUARY
16, 2021 ∂i∂t − αe + ( γ i + δ ) i − ∇ · ( nν i ∇ i ) = 0 (7) ∂r∂t − γ e e − γ i i − ∇ · ( nν r ∇ r ) = 0 (8) ∂d∂t − δi = 0 , (9)where β i and β e denote the transmission rates between symptomatic and susceptible individuals and asymptomatic andsusceptible individuals, respectively (units days − ), α denotes the incubation period (units days − ), γ e correspondsto the asymptomatic recovery rate (units days − ), γ i the symptomatic recovery rate (units days − ), δ represents themortality rate (units days − ), and ν s , ν e , ν i , ν r are the diffusion parameters of the different population groups asdenoted by the sub-scripted letters (units km persons − days − ). Note that all these parameters can be consideredtime and space-dependent.We may express the model (19)-(23) in the general form given by equation (1) by defining the tensors A , B , ν , and thevector f in the following way: A = α + γ e − α γ i + δ − γ e − γ i − δ (10) B = β e n s β i n s − β e n s − β i n s (11) ν = ν s ν e ν i ν r
00 0 0 0 0 (12) ν k = (cid:20) ν kxx ν kxy ν kyx ν kyy (cid:21) with k = s, e, i, r (13) f = . (14)Assuming additionally that the region of interest is isolated, we then prescribe the following homogeneous Neumannboundary conditions: ∇ s · n = 0 (15) ∇ e · n = 0 (16) ∇ i · n = 0 (17) ∇ r · n = 0 , (18)or simply ( ν · ∇ u ) · n = 0 . 4 PREPRINT - F
EBRUARY
16, 2021We note that one may also consider the following model related to (5)-(9): ∂s∂t + β i (cid:18) − An (cid:19) si + β e (cid:18) − An (cid:19) se − ∇ · ( nν s ∇ s ) = 0 (19) ∂e∂t − β i (cid:18) − An (cid:19) si − β e (cid:18) − An (cid:19) se + ( α + γ e ) e − ∇ · ( nν e ∇ e ) = 0 (20) ∂i∂t − αe + ( γ i + δ ) i − ∇ · ( nν i ∇ i ) = 0 (21) ∂r∂t − γ e e − γ i i − ∇ · ( nν r ∇ r ) = 0 (22) ∂d∂t − δi = 0 . (23)The model (19)-(23) is known as the density-dependent formulation and was studied in [8, 7]. The difference betweenthe models can be seen in equations (5)-(9) and (19)-(23). In (5)-(9), the contact terms are normalized by the livingpopulation, whereas such a normalization does not occur in (19)-(23). In the frequency-dependent formulation (5)-(9),this normalization implies that the contagion is independent of population density, while this is not the case in thedensity-dependent formulation, as the name may suggest [28, 29]. Both models may be valid and deliver accurateresults, depending on the physical situation, and we will show computations performed with both formulations in thepresent work.The other major difference between (5)-(9) and (19)-(23) is the presence of the Allee effect A . This term has been usedextensively in other settings, with the form used above inspired directly by applications in cancer modeling [30, 31].The Allee term serves to reduce transmission in areas where the population density is below a certain threshold A , bybringing the population in the exposed compartment to the susceptible compartment in such regions. Consequently, inthose areas, the population in compartments e and i tends to zero, which eventually cancels out the transfer term wherethe Allee term participates. In [8], this term served to accentuate the contagion in major urban regions while reducing itin less-populated areas, consistent with the observed physics. We briefly note also that as s , i , and e are all less than n by definition, we do not expect blowup of this term, even for very small n .For the numerical solution of (5)-(9), (19)-(23), we discretize in space using a Galerkin finite element formulation.The resulting systems of equations are stiff, leading us to employ implicit methods for time integration. We applythe second-order Backward Differentiation Formula (BDF2) in all cases, which offers second-order accuracy whileremaining unconditionally stable. We additionally make use of an adaptive mesh refinement and coarsening strategy(AMR/C), allowing us to resolve multiple scales. One may find more details about the adopted methods in [7, 27].All simulations have been performed using the libMesh , a C++ FEM open-source software library for parallel adaptivefinite element applications [32]. libMesh also interfaces with external solver packages like PETSc [33] and Trilinos[34]. It is an excellent tool for programming the finite element method and can be used for one-, two-, and three-dimensional steady and transient simulations. This library provides native support for AMR/C, thus providing a naturalenvironment for the present study. The main advantage of libMesh is the possibility of focusing on implementing thespecific features of the modeling without worrying about adaptivity and code parallelization. Consequently, the effort tobuild a high performance computing code tends to be minimized. In the following section, we demonstrate the robustness of the models (5)-(9), (19)-(23) by simulating the COVID-19outbreak for three different geographic regions: Lombardy (Italy), the U.S. state of Georgia, and the state of Rio deJaneiro (Brazil). These regions cover three different continents, and each exhibits very different geographical features,population characteristics, economic, social, cultural factors, and government responses. By achieving good results foreach case, we aim at showing that the modeling approach shown here is robust to general cases and does not requireextensive modification or specific parameter tuning to achieve reasonable results for a given case.The three cases shown were chosen for a number of reasons. Notably, they represent the authors’ home regions, so thereis an obvious interest from this point-of-view. Beyond this point, however, the regions also have specific characteristicsthat make them remarkable case studies. Italy is one of the hardest-hit regions in Europe, and Lombardy represented5
PREPRINT - F
EBRUARY
16, 2021the epicenter of the outbreak in Italy, as well as the first area in a Western country to exhibit community spread. Duringthe first wave of the epidemic in Spring 2020, Lombardy was disproportionately impacted compared to other regions inItaly, in which the epidemic was less severe. It remains the hardest-hit region in Italy, even after the second wave ofinfections has been less geographically restricted. As Lombardy was studied specifically in [8, 7], it also represents anatural choice for validating the new code and the adaptive-meshing strategy implemented in the present work.The U.S. state of Georgia is a logical choice for several reasons. The state has 159 counties, second to only Texasamong U.S. states. However, the average county size for Georgia is less than half of Texas, giving the state a goodcounty-by-county spatial resolution, and allowing us to evaluate the spatial accuracy of the model in reasonable detail.Further, the state was the first among U.S. states to reopen in late April, and has since taken a large number of policymeasures, in both relaxing and increasing restrictions, providing us a large number of decisions to incorporate andanalyze in our model.Rio de Janeiro is also a sensible choice for many of the same reasons as Georgia. At the time of writing, Brazil is thethird-most severely impacted country in the world in terms of total cases, and the second one in terms of total deaths.Rio de Janeiro has the largest number of deaths per inhabitant within Brazil, making its simulation worthwhile onseverity alone. Aside from this aspect, the area also exhibits interesting geographical features, notably its naturalharbor and large population centers on either side. Additionally, there is good availability of mobility data, whichwe have incorporated into the diffusion parameters. The use of such data is logical for a model of this type, and itsconsideration represents an interesting and important component for this modeling approach.
In our first numerical test, we seek to reproduce, and possibly improve, the simulations of Lombardy, Italy, as shown in[8, 7]. This model was already shown in [8] to predict the first wave of the COVID-19 outbreak in Lombardy with goodaccuracy. We note that for this simulation we use the density-dependent model with the Allee term (19)-(23). Whileusing the same parameters as shown in [8, 7], we incorporate two major improvements. First, we employ the adaptivemesh refinement and coarsening strategy shown in [27]. Second, rather than using Backward-Euler time-stepping, weutilize a BDF2 scheme, as [35, 7] suggest that such a scheme is better suited to the system of equations in question.Thus, the primary motivation of this simulation is the validation of the current code and algorithmic framework.The AMR/C procedure uses a local error estimator to drive the refinement and coarsening procedure, considering theerror of an element relative to its neighbor elements in the mesh. The elements are flagged depending on a refining( r f ) and a coarsening ( c f ) fraction. The refinement level is limited by a maximum h -level ( h max ). In this case, for theAMR/C procedure, we set h max = 1 , r f = 0 . , c f = 0 . . We apply the adaptive mesh refinement every 4 time-steps.The original mesh has 10,987 linear triangular elements, and, after refinement, the minimum spatial resolution is about1 kilometer. We initially refine the domain in one level, and the time-step is defined as ∆ t = 0 . days .In Fig. 1 we compare the results obtained here with those from [7]. We observe the same qualitative dynamics; theepidemic follows the same path, moving from the southern areas of the region northward into the major metropolitanareas, and the Milan area in particular. The changes to the meshing and time-stepping schemes result in somewhat lowercontagion, particularly in the areas in and around the city of Milan in the western part of the region. Examining theprogression of the mesh adaptation in Fig. 1, we see that the AMR/C algorithm successfully follows the transmissionpath, changing in time to adapt to the evolving geography of the epidemic in time. We briefly note that when usingthe current code with Backward-Euler and a fixed mesh (not shown), the results were identical to those obtained in[8]; thus, we can conclude that the differences between the solution depicted in the present work and that shown in [8]represent improvements resulting from the differences in meshing and time-stepping used in the present one. From this,we conclude that the current code and algorithmic framework is validated. From now on we use the frequency-dependent model (5)-(9) together with adaptive mesh refinement and the BDF2time-stepping scheme. In the next test, we aim at reproducing the COVID-19 outbreak in the U.S. state of Georgia .
We have first obtained the map of the state of Georgia (GA) along with the county boundaries in shapefile formatfrom https://arc-garc.opendata.arcgis.com/datasets/dc20713282734a73abe990995de40497_68 . Totriangulate the GA region, we follow the steps given by [14]:• Load the GA map file in the freely available QGIS software [36] ;6
PREPRINT - F
EBRUARY
16, 2021Figure 1: Spread of the COVID-19 virus spread at Lombardy, Italy. Comparison between the results obtained in [8], inwhich the time integration was done using the Backward-Euler method with a fixed mesh, and the present work, inwhich we use the BDF2 method and AMR.• Coarse grain the outer boundary segments using the Simplify tool in QGIS. The original map has some verysmall length segments, which may create problems in triangulation or result in a very fine mesh;• Obtain the vertices using the Extract Vertices tool in QGIS and save the vertices layer using the save layer asoption. Select As XY in the Graphical category and save the file in a .csv format;• Prepare a Gmsh input file using the vertices’ file for triangulation.In Figure 2, we show the generated grid for the state of Georgia. The original mesh has 32,056 linear triangularelements, and after refinement, the minimum spatial resolution is about 2 kilometers. We initially refine the domain inone level. For the AMR/C procedure, we set h max = 1 , r f = 0 . , c f = 0 . . We apply the adaptive mesh refinementevery 4 time-steps. The time-step is defined as ∆ t = 0 . days .We define the initial infected population accordingly to the data provided by the Johns Hopkins University [37] . Wedefine the beginning of the simulation on 25 March 2020 and simulate 250 days. The exposed population is moredifficult to estimate, as they consist of asymptomatic and pre-symptomatic individuals. Here, we consider 10 times thenumber of infected [38]. The susceptible population is based on the estimation of the population of each county, givenby the U.S. Census Bureau[39]. The populations are divided by the area of each county and distributed on the 159areas of the mesh as people/km . In Figs. 3, 4, and 5 we show the initial conditions.The parameters of the simulation are defined as: α = 1 / day − , γ i = 1 / day − , γ e = 1 / day − , δ = 1 / day − . These values are based on available data from the literature regarding the mortality, incubation period, andrecovery time for infected and asymptomatic patients [8].The contact rate and diffusion are estimated based on the policies adopted in the Georgia state during the period ofthe simulation ( https://coronavirus.jhu.edu/data/state-timeline/new-deaths/georgia/ ). This leadsto the time-dependent values as shown in Fig. 6. 7 PREPRINT - F
EBRUARY
16, 2021Figure 2: Map of the state of GA partitioned into 159 internal counties.Figure 3: Initial susceptible population ( people/km ).8 PREPRINT - F
EBRUARY
16, 2021Figure 4: Initial exposed population ( people/km ).Figure 5: Initial infected population ( people/km ).9 PREPRINT - F
EBRUARY
16, 2021Figure 6: Diffusion (left) and contact rate (right) in time for Georgia.
To compare our model results with the data available, we compare the number of deaths caused by the COVID-19. First,we show the values of the whole state. Then, we integrate the values of some counties and plot them in time (Fig. 8 and9). Fig. 10 also shows some snapshots of the simulation with the populations and mesh at different time-steps.Referring to the entire state in Fig. 7, we observe that the model accurately captures both qualitative and quantitativedynamics of the epidemic with good accuracy. The difference in relative L norm between the measured and simulateddeaths is 5.75%, confirming the agreement. We note that a small discrepancy between the simulated and measured dataoccurs at the end of the simulated period; this is caused by a series of deaths being added retroactively to the data, withthe actual date of these deaths being unknown.Turning our attention to the individual counties, we see more variation in the results. We start by noting that the initialdynamics (first 60 days) is predicted well nearly everywhere, with a divergence between the simulated and measureddeaths becoming more apparent as the epidemic progresses. This divergence generally indicates that such models’predictive power may decrease in time, as one may reasonably expect. Nonetheless, we do observe encouraging resultsover the entire simulated time frame, particularly in certain areas. We first assess counties comprising the Atlantametropolitan area. The heavily populated counties of Dekalb, Gwinett, Clayton, and Cobb show very good agreementwith the simulated data. The most populous county in Georgia, Fulton County, is over-predicted. On the fringe of theAtlanta metropolitan area, Forsyth, Clarke, and Cherokee counties were also over-predicted.Outside of the Atlanta area, we found a general underestimation of the outbreak severity in the Bibb, Richmond, andChatham counties. As the “first wave" of infections is well-simulated, this likely represents a shortcoming of thediffusive model. In particular, after the initial re-openings, additional source terms may need to be included in certainareas to model the effects of travelers arriving from elsewhere, as our model does not account for nonlocal movements.In contrast, Dougherty County, which was heavily impacted early in the epidemic, shows the opposite trend, and theoutbreak during the later stages of the simulated period is less severe than as predicted by the model. This shows thatthe model is sensitive to the initial conditions, and spatially-variable parameters, in particular contact rates, may beneeded to account for such phenomena properly. The last simulation aims at reproducing the COVID-19 outbreak in the state of Rio de Janeiro, Brazil. Rio de Janeirohas the second number of deaths and infected in Brazil, but the largest number of deaths/inhabitant as per January2021 We have obtained the map of the state of Rio de Janeiro (RJ) along with the county boundaries inshapefile format from ftp://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/ PREPRINT - F
EBRUARY
16, 2021Figure 7: Comparison between simulation and real data of deaths at Georgia (total). malhas_municipais/municipio_2017/UFs/ . To triangulate the RJ region, we follow the same steps explained inthe previous section.In Fig. 11, we show the grid generated for the state of Rio de Janeiro. We eliminated some islands to facilitatethe simulation. The original mesh has 11,632 linear triangular elements, and after refinement, the minimum spatialresolution is about 500 meters. We initially refine the domain in one level. For the AMR/C procedure, we set h max = 1 , r f = 0 . , c f = 0 . . We apply the adaptive mesh refinement every 4 time-steps. The time-step is defined as ∆ t = 0 . days .We define the initial infected population accordingly to the data provided by the Brazilian Ministry of Health, collectedby https://covid19br.wcota.me/ . We define the beginning of the simulation on 25 March 2020 and simulate 180days. However, since there is a delay until people develop symptoms, another delay in receiving the test results, andanother one for the case being disclosed, we consider that the initial infected population is equal to the case numbersprovided for each county on 1 April 2020. As in the case of Georgia, we initialize the exposed population as 10 timesthe number of infected [38]. The susceptible population is based on the estimation of the population of each county,given by the Brazilian Institute for Geography and Statistics . The populations are divided by the area of each countyand distributed on the 92 areas of the mesh as people/km .In Figs. 12, 13, and 14 we show the initial conditions. The population of the state of Rio de Janeiro is concentrated nearthe capital and the metropolitan region.The parameters of the simulation are defined as: α = 1 / day − , γ i = 1 / day − , γ e = 1 / day − , δ = 1 / day − . These values are based on available data from the literature regarding the mortality, incubation period, andrecovery time for infected and asymptomatic patients [8].The contact rate is estimated based on the reproduction number estimation given in https://perone.github.io/covid19analysis/ . We define β e = β i = 0 . days − , and we multiply this value by the time functiongiven by Fig. 15b, which is a simplification of the behavior of the reproduction number between 25 March 2020and 21 September 2020. The diffusion coefficient is based on the social distancing estimation given by https://coronavirus.ufrj.br/covidimetro/ , which presents the perception of average weekly social isolation basedon the movement of cell phones in the state of RJ. The diffusion behaves opposite to the social isolation, i.e.,when one increases, the other one decreases. Therefore, we set ν s = ν e = ν r = 1 × − and ν i = 1 × − km persons − days − , and scale these values by a time-dependent function (Fig. 15a) that tries to represent thesocial isolation during the 180 days of simulation. PREPRINT - F
EBRUARY
16, 2021Figure 8: Comparison between simulation and real data of deaths at Georgia (per county) - part 1.12
PREPRINT - F
EBRUARY
16, 2021Figure 9: Comparison between simulation and real data of deaths at Georgia (per county) - part 2.13
PREPRINT - F
EBRUARY
16, 2021Figure 10: Snapshots of Georgia simulation. Top to bottom: Susceptible, Exposed, Infected, Recovery and Deceasedpopulation ( people/km ), followed by the mesh at different time-steps.14 PREPRINT - F
EBRUARY
16, 2021Figure 11: Map of the state of RJ partitioned into 92 internal counties.
To compare the results of our model with the data available, we compare the number of deaths caused by the COVID-19.First, we show the values of the whole state (Fig. 16). Then, we integrate the values of some counties and plot them intime (Figs. 17 and 18). Fig. 19 also shows some snapshots of the simulation with the populations and mesh at differenttime-steps.As in the previous case in Georgia, most selected counties show good agreement with the available data. Notably, overthe entire region, the relative difference in L norm between the observed and simulated fatalities is only 10.06%. Overthe first several weeks, the results are in good agreement nearly everywhere, and, in the areas around the city of Rio deJaneiro, we capture the dynamics very accurately over the entire considered time interval. In addition to the countyof Rio de Janeiro, we observe particularly good agreement in the highly populated areas of São Gonçalo, Duque deCaxias, Belford Roxo, and Itaborai. In the areas of Petropolis, Volta Redonda, and Nova Igauçu, we obtain somewhatless accurate, but still nonetheless reasonable results. In these latter cases, we observe good agreement over the initialphases of the dynamics, with some differences later in the time interval. These deviations are unsurprising, as one maynaturally expect increased difficulty further in time.We do observe some notable instances of overestimation in Niterói and underestimation in Cabo Frio and Camposdos Goytacazes. In the latter cases, there were no initial infected or exposed population in this simulation, makingdiffusion the only possible way for the virus to reach the areas. Given the large distances between these areas andthe hotspots near the city center of Rio de Janeiro (157 km to Cabo Frio and 279 km to Campos dos Goytacazes), aswell as large sparsely-populated areas between them, population-weighted diffusion was not able to properly representthese dynamics. The model best captures spatial dynamics over relatively short distances overpopulated regions. Thisrepresents a shortcoming of the model, and additional terms, such as source/sink terms, fractional diffusion operators,or bilaplacian diffusion terms, may be required to properly account for such nonlocal dynamics. In the case ofNiterói, the overprediction may result from an overestimation of the initial exposed population or failure to take intoaccount specific local policies. This is similar to what we observed in Georgia, in which certain particular regions had15 PREPRINT - F
EBRUARY
16, 2021Figure 12: Initial susceptible population ( people/km ).Figure 13: Initial exposed population ( people/km ).16 PREPRINT - F
EBRUARY
16, 2021Figure 14: Initial infected population ( people/km ). (a) Function for ν . (b) Function for β . Figure 15: Functions that multiplies the diffusion (left) and contact rate (right) in time for Rio de Janeiro.17
PREPRINT - F
EBRUARY
16, 2021Figure 16: Comparison between simulation and real data of deaths at Rio de Janeiro state (total).dynamics that deviated from the measured data. However, as in Georgia’s case, we observe very good agreement whenconsidering the totality of the region, with particular areas less well-represented due to the natural limitations of ourmodeling approach. In these instances, such cases provide clear directions in which we may improve the modelingframework.
In this work, we have established the robustness of the SEIRD model introduced in [8] and further extended in [7, 27, 14]by examining both its frequency- and density-dependent formulations and applying it to three settings over differentcontinents: the region of Lombardy, Italy (density-dependent), the U.S. state of Georgia, and the state of Rio de Janeiro,Brazil (frequency-dependent). All these regions have very different geographic characteristics, population densitypatterns, policy response measures, and cultural contexts generally; despite this, the models were able to adequatelyreproduce the spatio-temporal contagion dynamics with minimal parameter differences and very little tuning. Thisprovides strong evidence towards the robustness of the model, as it can produce adequate results across different settingswithout requiring extensive parameter fitting and learning.The results obtained showed good qualitative agreement generally across all regions; however, they show clear room forimprovement. In the cases of Georgia and Rio de Janeiro, we find that the initial conditions have a large influence,even over the long-term dynamics of the system. This is particularly apparent in the case of Dougherty County inGeorgia, where an early outbreak led to a long-term overestimation of contagion, and in Campos dos Goytacazes andCabo Frio in Rio de Janeiro, where a relative lack of early cases led to an under-prediction of the long term contagioneffects. This is likely due to the diffusion formulation limitations, which cannot correctly account for non-local mobilityacross different areas, such as returning travelers. Such effects could be incorporated in a number of ways by includingmore sophisticated source/sink terms, fractional diffusion operators, or bilaplacian terms. In Lombardy, we observean under-prediction of contagion in less-populated regions; this appears to be characteristic of the density-dependentformulation, and the frequency-dependent models applied in Georgia and Rio de Janeiro do not suffer from this problem.Across all cases, the same parameter values were used everywhere. While this represents evidence towards robustnessand is in some sense a positive aspect, results could, of course, be further improved by using different definitionsin different regions, something explored briefly in [8] and in more detail in [14]. Such approaches may incorporatemore sophisticated methods, including optimization and machine-learning techniques. These shortcomings representnatural directions for future work in this area. Nonetheless, we believe that our results clearly establish the viability androbustness of this PDE modeling framework, suggesting this as a suitable path for future research, and ultimately mayperhaps prove useful to public health decision-makers. 18
PREPRINT - F
EBRUARY
16, 2021Figure 17: Comparison between simulation and real data of deaths at RJ (per county) - part 1.19
PREPRINT - F
EBRUARY
16, 2021Figure 18: Comparison between simulation and real data of deaths at RJ (per county) - part 2.20
PREPRINT - F
EBRUARY
16, 2021Figure 19: Snapshots of Rio de Janeiro simulation. Top to bottom: Susceptible, Exposed, Infected, Recovery andDeceased population ( people/km ), followed by the mesh at different time-steps.21 PREPRINT - F
EBRUARY
16, 2021
This research was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil(CAPES) - Finance Code 001 and CAPES TecnoDigital Project 223038.014313/2020-19. This research has alsoreceived funding from CNPq and FAPERJ. A. Reali was partially supported by the Italian Ministry of University andResearch (MIUR) through the PRIN project XFAST-SIMS (No. 20173C478N).
References [1] Richard A Erickson, Steven M Presley, Linda JS Allen, Kevin R Long, and Stephen B Cox. A dengue model witha dynamic aedes albopictus vector population.
Ecological Modelling , 221(24):2899–2908, 2010.[2] Eber Dantas, Michel Tosin, and Americo Cunha Jr. Calibration of a seir–sei epidemic model to describe the zikavirus outbreak in brazil.
Applied Mathematics and Computation , 338:249–259, 2018.[3] Zindoga Mukandavire, Prasenjit Das, Christinah Chiyaka, and Farai Nyabadza. Global analysis of an hiv/aidsepidemic model.
World Journal of Modelling and Simulation , 6(3):231–240, 2010.[4] Chris Dye and Nigel Gay. Modeling the sars epidemic.
Science , 300(5627):1884–1885, 2003.[5] Alemayehu Midekisa, Gabriel Senay, Geoffrey M Henebry, Paulos Semuniguse, and Michael C Wimberly.Remote sensing-based time series models for malaria early warning in the highlands of ethiopia.
Malaria Journal ,11(1):165, 2012.[6] Phenyo E Lekone and Bärbel F Finkenstädt. Statistical inference in a stochastic epidemic seir model with controlintervention: Ebola as a case study.
Biometrics , 62(4):1170–1177, 2006.[7] Alex Viguerie, Alessandro Veneziani, Guillermo Lorenzo, Davide Baroli, Nicole Aretz-Nellesen, Alessia Patton,Thomas E Yankeelov, Alessandro Reali, Thomas JR Hughes, and Ferdinando Auricchio. Diffusion–reactioncompartmental models formulated in a continuum mechanics framework: application to covid-19, mathematicalanalysis, and numerical study.
Computational Mechanics , pages 1–22, 2020.[8] Alex Viguerie, Guillermo Lorenzo, Ferdinando Auricchio, Davide Baroli, Thomas JR Hughes, Alessia Patton,Alessandro Reali, Thomas E Yankeelov, and Alessandro Veneziani. Simulating the spread of covid-19 viaspatially-resolved susceptible-exposed-infected-recovered-deceased (seird) model with heterogeneous diffusion.
Applied Mathematics Letters , 111:106617, 2021.[9] Julien Arino and Stéphanie Portet. A simple model for covid-19.
Infectious Disease Modelling , 5:309–315, 2020.[10] Giulia Giordano, Franco Blanchini, Raffaele Bruno, Patrizio Colaneri, Alessandro Di Filippo, Angela Di Matteo,and Marta Colaneri. Modelling the covid-19 epidemic and implementation of population-wide interventions initaly.
Nature Medicine , pages 1–6, 2020.[11] José M Carcione, Juan E Santos, Claudio Bagaini, and Jing Ba. A simulation of a covid-19 epidemic based on adeterministic seir model.
Frontiers in Public Health , 8:230, 2020.[12] Diego Tavares Volpatto, Anna Claudia Mello Resende, Lucas Anjos, Joao Vitor Oliveira Silva, Claudia MazzaDias, Regina Cerqueira Almeida, and Sandra Mara Cardoso Malta. Spreading of covid-19 in brazil: Impacts anduncertainties in social distancing strategies. medRxiv , 2020.[13] Juliane F Oliveira, Daniel CP Jorge, Rafael V Veiga, Moreno S Rodrigues, Matheus F Torquato, Nivea B da Silva,Rosemeire L Fiaccone, Luciana L Cardim, Felipe AC Pereira, Caio P de Castro, et al. Mathematical modeling ofcovid-19 in 14.8 million individuals in bahia, brazil.
Nature communications , 12(1):1–13, 2021.[14] Prashant K Jha, Lianghao Cao, and J Tinsley Oden. Bayesian-based predictions of covid-19 evolution in texasusing multispecies mixture-theoretic continuum models.
Computational Mechanics , 66(5):1055–1068, 2020.[15] Ivan Korolev. Identification and estimation of the seird epidemic model for covid-19.
Journal of econometrics ,220(1):63–85, 2020.[16] Fleurianne Bertrand and Emilie Pirch. Least-squares finite element method for a meso-scale model of the spreadof covid-19.
Computation , 9(2), 2021.[17] Fred Brauer, Carlos Castillo-Chavez, and Zhilan Feng.
Mathematical models in epidemiology . Springer, 2019.[18] Marino Gatto, Enrico Bertuzzo, Lorenzo Mari, Stefano Miccoli, Luca Carraro, Renato Casagrandi, and AndreaRinaldo. Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures.
Proceedings of the National Academy of Sciences , 117(19):10484–10491, 2020.22
PREPRINT - F
EBRUARY
16, 2021[19] Francesc Aràndiga, Antonio Baeza, Isabel Cordero-Carrión, Rosa Donat, M Carmen Martí, Pep Mulet, andDionisio F Yáñez. A spatial-temporal model for the evolution of the covid-19 pandemic in spain includingmobility.
Mathematics , 8(10):1677, 2020.[20] Valerio La Gatta, Vincenzo Moscato, Marco Postiglione, and Giancarlo Sperli. An epidemiological neural networkexploiting dynamic graph structured data applied to the covid-19 outbreak.
IEEE Transactions on Big Data , 2020.[21] TI Zohdi. An agent-based computational framework for simulation of global pandemic and social response onplanet x.
Computational Mechanics , 66(5):1195–1209, 2020.[22] Petrônio CL Silva, Paulo VC Batista, Hélder S Lima, Marcos A Alves, Frederico G Guimarães, and Rodrigo CPSilva. Covid-abs: An agent-based model of covid-19 epidemic to simulate health and economic effects of socialdistancing interventions.
Chaos, Solitons & Fractals , 139:110088, 2020.[23] Navid Mahdizadeh Gharakhanlou and Navid Hooshangi. Spatio-temporal simulation of the novel coronavirus(covid-19) outbreak using the agent-based modeling approach (case study: Urmia, iran).
Informatics in MedicineUnlocked , 20:100403, 2020.[24] Joshua P Keller, Luca Gerardo-Giorda, and Alessandro Veneziani. Numerical simulation of a susceptible–exposed–infectious space-continuous model for the spread of rabies in raccoons across a realistic landscape.
Journal ofBiological Dynamics , 7(sup1):31–46, 2013.[25] Mi-Young Kim. Galerkin methods for a model of population dynamics with nonlinear diffusion.
NumericalMethods for Partial Differential Equations: An International Journal , 12(1):59–73, 1996.[26] Robert Stephen Cantrell and Chris Cosner.
Spatial ecology via reaction-diffusion equations . John Wiley & Sons,2004.[27] Malú Grave and Alvaro L. G. A. Coutinho. Adaptive mesh refinement and coarsening for diffusion-reactionepidemiological models, 2020.[28] Madan K Oli, Meenakshi Venkataraman, Paul A Klein, Lori D Wendland, and Mary B Brown. Populationdynamics of infectious diseases: a discrete time model.
Ecological Modelling , 198(1-2):183–194, 2006.[29] Peter H Thrall, Arjen Biere, and Marcy K Uyenoyama. Frequency-dependent disease transmission and thedynamics of the silene-ustilago host-pathogen system.
The American Naturalist , 145(1):43–62, 1995.[30] Kaitlyn E Johnson, Grant Howard, William Mo, Michael K Strasser, Ernesto ABF Lima, Sui Huang, and AmyBrock. Cancer cell population growth kinetics at low densities deviate from the exponential growth model andsuggest an allee effect.
PLoS biology , 17(8):e3000399, 2019.[31] Marcello Delitala and Mario Ferraro. Is the Alee effect relevant in cancer evolution and therapy?
AIMSMathematics , 5(6):7649–7660, 2020.[32] B. S. Kirk, J. W. Peterson, R. H. Stogner, and G. F. Carey. libmesh: a c++ library for parallel adaptive meshrefinement/coarsening simulations.
Journal Engineering with Computers , 22(3):237–254, 2006.[33] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin,Alp Dener, Victor Eijkhout, William D. Gropp, Dmitry Karpeyev, Dinesh Kaushik, Matthew G. Knepley, Dave A.May, Lois Curfman McInnes, Richard Tran Mills, Todd Munson, Karl Rupp, Patrick Sanan, Barry F. Smith,Stefano Zampini, Hong Zhang, and Hong Zhang. PETSc Web page. , 2019.[34] The Trilinos Project Team. The Trilinos Project Website.[35] Malú Grave, José J. Camata, and Alvaro L. G. A. Coutinho. A new convected level-set method for gas bubbledynamics.
Computers & Fluids , 209:104667, 2020.[36] QGIS Development Team.
QGIS Geographic Information System . Open Source Geospatial Foundation, 2009.[37] Ensheng Dong, Hongru Du, and Lauren Gardner. An interactive web-based dashboard to track covid-19 in realtime.
The Lancet infectious diseases , 20(5):533–534, 2020.[38] Ruy Freitas Reis, Bárbara de Melo Quintela, Joventino de Oliveira Campos, Johnny Moreira Gomes, Bernardo Mar-tins Rocha, Marcelo Lobosco, and Rodrigo Weber dos Santos. Characterization of the covid-19 pandemic and theimpact of uncertainties, mitigation strategies, and underreporting of cases in south korea, italy, and brazil.
Chaos,Solitons & Fractals , page 109888, 2020.[39] United States Census Bureau.