key: cord-0905712-nbdc9kkz authors: Viguerie, Alex; Lorenzo, Guillermo; Auricchio, Ferdinando; Baroli, Davide; Hughes, Thomas J.R.; Patton, Alessia; Reali, Alessandro; Yankeelov, Thomas E.; Veneziani, Alessandro title: Simulating the spread of COVID-19 via a spatially-resolved susceptible–exposed–infected–recovered–deceased (SEIRD) model with heterogeneous diffusion date: 2020-07-15 journal: Appl Math Lett DOI: 10.1016/j.aml.2020.106617 sha: 893d33bc11f776323afb3db92f2751895a4497fe doc_id: 905712 cord_uid: nbdc9kkz We present an early version of a Susceptible–Exposed–Infected–Recovered–Deceased (SEIRD) mathematical model based on partial differential equations coupled with a heterogeneous diffusion model. The model describes the spatio-temporal spread of the COVID-19 pandemic, and aims to capture dynamics also based on human habits and geographical features. To test the model, we compare the outputs generated by a finite-element solver with measured data over the Italian region of Lombardy, which has been heavily impacted by this crisis between February and April 2020. Our results show a strong qualitative agreement between the simulated forecast of the spatio-temporal COVID-19 spread in Lombardy and epidemiological data collected at the municipality level. Additional simulations exploring alternative scenarios for the relaxation of lockdown restrictions suggest that reopening strategies should account for local population densities and the specific dynamics of the contagion. Thus, we argue that data-driven simulations of our model could ultimately inform health authorities to design effective pandemic-arresting measures and anticipate the geographical allocation of crucial medical resources. The outbreak of COVID-19 in 2020 has caused widespread disruption throughout the world, leading to substantial damage in terms of both human lives and economic cost. To arrest the spread of the disease, governments have enacted unprecedented measures, including quarantines, curfews, lockdowns, and suspension of travel. The wide-reaching ramifications of such measures, deemed by many experts as necessary, are driven in part by a lack of clear information about the spatio-temporal spread of COVID-19. Indeed, the absence of reliable data regarding disease transmission has necessarily led to cautious responses. These recent events, which have required important decisions based on forecasts, have demonstrated more than ever the need for reliable tools intended to model the spread of COVID-19 and other infectious diseases [17] . A particularly urgent need is the geo-localization of outbreaks, as this may allow a more effective allocation of medical resources. Several notable models of this outbreak have been presented; indeed, at the time of this writing there are over 1,000 COVID-19 articles on MedRXiv, many of which address the modeling of disease spread. Some models aim at offering specific evaluations of policy responses based on the implementation of different social distancing measures with combined compartmental and empirical approaches [5, 7, 8] . Rather than adopting a deterministic, mechanism-based model, Zhang et al. employed a statistical approach to analyze the spatio-temporal dynamics of COVID-19 [20] . In Gatto et al. a combined statistical and compartmental approach was employed, in which spatial dependence is addressed by dividing the region of interest (Italy) into local communities connected by a network structure [7] . Here, we propose an alternative approach, using a partial-differential-equation (PDE) model designed to capture the continuous spatio-temporal dynamics of COVID-19. We leverage a compartmental SEIRD (susceptible, exposed, infected, recovered, deceased ) model that incorporates the spatial spread of the disease with inhomogeneous diffusion terms [9, 10, 12, 13] . The rationale is that the diffusion operator, properly tuned to account for local natural or social inhomogeneities (e.g., mountains, rivers, highways) may describe the local movement of the different populations in a deterministic way, as the limit of a Brownian motion [18] . This is critical to accurately account for information relevant to the outbreak dynamics, such as local population densities, which vary in space and time. While a mathematical description of non-local dynamics is still possible in terms of fractional differential operators [11] , we postpone this approach to a follow-up of the present work. Hence, our modeling approach is more appropriate for the local dynamics on mesoscales, such as regions within Italy. Thus, to evaluate the model efficacy, we run a simulation study of the COVID-19 outbreak in the Italian region of Lombardy, which has been severely impacted by the COVID-19 crisis between February and April of 2020 and for which the necessary data was available. Also, the high density of Lombardy's population and transportation network is specifically suitable to our modeling approach. Our simulations show a remarkable qualitative agreement with the reported epidemiological data. We further explore various reopening scenarios, obtaining contrasting results that highlight the importance of considering local population densities and contagion dynamics. The paper outline is as follows. In Section 2, we describe the SEIRD model. Then, Section 3 addresses the numerical implementation of our model and Section 4 presents the results of the simulation study in Lombardy. We conclude in Section 5 by examining the shortcomings observed from our simulations and discussing the additional work required to improve model accuracy and practical relevance. Let Ω ⊂ R 2 be a simply connected domain of interest and [0, T ] a generic time interval. We denote the densities of the susceptible, exposed, infected, recovered and deceased populations as s(x, t), e(x, t), i(x, t), r(x, t), and d(x, t) respectively. Also, let n(x, t) denote the sum of the living population; i.e., n(x, t) = s(x, t) + e(x, t) + i(x, t) + r(x, t). Then, our model is comprised of the following system of coupled PDEs over Ω × [0, T ]: where α is the birth rate, σ is the inverse of the incubation period, φ e is the asymptomatic recovery rate, φ r is the infected recovery rate, φ d is the infected mortality rate, β e is the asymptomatic contact rate, β i is the symptomatic contact rate, µ is the general (non-COVID-19) mortality rate, and ν s , ν e , ν i , and ν r are diffusion parameters respectively corresponding to the different population groups. Each of these parameters may depend on time, space, or the model compartments. We also consider the Allee effect (depensation), characterized by the parameter A. In this particular setting, the Allee effect serves to model the tendency of outbreaks to cluster towards large population centers. Specific parameter selection as well as initial and boundary conditions are discussed in Section 3. Fig. 1 shows the dynamics of contagion between the compartments in our model. We remark that our model accounts for asymptomatic transmission, which is considered a pivotal driver of the COVID-19 pandemic [3, 8, 15] . Eqs. (1)-(2) show that exposed asymptomatic patients may transmit COVID-19 to susceptible individuals at contact rate β e . This aligns with recent studies suggesting that patients may transmit COVID-19 almost immediately after exposure [3, 8, 15] . Additionally, Eqns. (2) and (4) involve a fraction φ e of exposed patients that do not develop symptoms and move directly into the J o u r n a l P r e -p r o o f Journal Pre-proof Figure 1 : Flow chart describing the dynamics of contagion between the population subgroups considered in our model recovered population. We also assume that recovered patients are immune, as we do not include any backflow from Eq. (4) to Eq. (1). (We note that this is a current source of debate, but consistent with the existing literature for the time scale of months considered here [14] ). The spatial movement over a large population is described by an inhomogeneous random walk, which in the limit tends to a second order differential operator [18] . The diffusivity coefficient is proportional to the population and can be locally adjusted to incorporate geographical or human-related inhomogeneities [10] . We use a finite-element spatial discretization of the Italian region of Lombardy, consisting of an unstructured mesh containing 30,407 triangles (a mesh convergence analysis was first performed; data not shown). We use the backward-Euler method for time integration and solve each time step fully implicitly with a Picard iteration for stability. The resulting linear systems are solved by the GMRES algorithm using a Jacobi preconditioner. Initial conditions for the subpopulations s, e, i, r, and d in the model are defined by means of Gaussian circular functions centered at the latitude and longitudinal coordinates of each municipality with 10,000 or more inhabitants, and weighted by the municipality's population size and geographic area. These initial conditions correspond to the data provided by Lab24 [1] from the date 27 February 2020, featuring a severe outbreak in the province of Lodi, and moderate numbers of exposed and infected individuals in the provinces of Bergamo, Brescia, and Cremona (see Fig. 2 ). We used homogeneous Neumann boundary conditions for simplicity, mimicking a complete isolation of the region. We acknowledge that this is likely unrealistic (e.g. if lockdown orders are relaxed) and will be considered in future efforts. We assume σ=1/7 day −1 , φ r =1/24 day −1 , φ d =1/160 day −1 , and φ e =1/6 day −1 . These values were based on available data from the literature regarding the mortality, incubation period, and recovery time for infected and asymptomatic patients [2, 3, 5, 8, 15] . Additionally, we do not consider births or non-COVID-19 mortality (i.e., we set α = 0 and µ = 0, respectively), given the time scale of months in our simulations. The remainder of the model parameters are estimated in a two-step approach. First, we fit a 0D SEIRD version of our model (i.e., consisting of a system of exclusively time-dependent ordinary differential equations and no diffusion terms) to match the temporal dynamics of the outbreak. Then, we iteratively refined these values by means of recursive simulations using Eqs. (1)-(5) to match the spatiotemporal epidemiological data. We use the R 2 coefficient and the root mean squared error (RMSE) to assess the goodness-of-fit. Given the uncertainty in the currently available COVID-19 data, we think that parameter estimation aiming at matching the dynamics of all model compartments is not viable. As not every member of the population is tested for infection and asymptomatic cases are known to exist in possibly large numbers, we think that the available data of infected cases might lead to unrealistic parameter fitting. Conversely, the data reported for COVID-19 deaths offer more reliability to calibrate the model parameters. Therefore, we pursue quantitative agreement in the deceased compartment (i.e., d), and qualitative agreement for the rest of the model subgroups (i.e., s,e,i,r). In Section 4, our results will focus on the model forecasts of exposed and infected cases because these are key data for public health officials, e.g., in deciding resource allocation and measures to prevent contagion. 0.0198, 0.0090, and 0.0075 km 2 · persons −1 · day −1 over the respective phases. The Allee term A is set to 1,000 persons, and we fix ν i = 1.0 · 10 −4 km 2 · persons −1 · day −1 throughout, assuming that symptomatic individuals are largely immobile. Lombardy Fig. 2 shows the evolving spatial pattern of the COVID-19 outbreak in Lombardy, beginning with exposure in Bergamo, Brescia, Cremona and Lodi. The contagion moves north from Lodi into Milan via the southern suburbs and eventually reaches the city center. We note that although Lodi and Cremona are the most affected areas at the onset of the outbreak, they quickly improve and avoid the explosive growth found in Milan, Bergamo and Brescia. This is consistent with the reported data [1] . However, we found Cremona to be somewhat underpredicted by our model. This might be attributable to the presence of the neighboring city of Piacenza, which shares its metropolitan area with Cremona but is not included in our simulations because it belongs to the adjacent region of Emilia-Romagna. The inclusion of a possible inflow of infected through the southern border with Emilia-Romagna does show a small spike in Cremona, as reported in Fig. 3 . In this case, in the absence of available data, we postulated with an educated guess an abnormal influx from the South (100 susceptible persons Km −1 day −1 ), 1 exposed person Km −1 day −1 , 1 infected person Km −1 day −1 ). This result pinpoints the importance of an accurate calibration of the boundary conditions. In Fig. 2 , we demonstrate the remarkable qualitative agreement in the outbreak dynamics between our model forecasts and data in the three main affected areas: Milan, Bergamo, and Brescia. The R 2 between the model forecast and the data of infected cases are 0.997, 0.977, 0.976, and 0.998 for all Lombardy, Bergamo, Brescia, and Milan, respectively. We observe that the outbreak emerges in Milan later, where it grows more steadily, eventually becoming the most affected area in Lombardy. We also note that the lockdowns appear to have effectively halted the spread in Bergamo and Brescia. These restrictions notably reduce the spread in Milan, limiting the virus to a linear growth pattern, but fail stop it. We observe that our simulations predict a larger number of infections than the reported data. This results from using the data for deceased cases for calibration, which is comparatively more accurate than infections (see Section 3). We obtained R 2 = 0.972 and range-normalized RMSE=7.6% for this subgroup. Thus, the difference in predicted and measured infections suggests a lack in the reporting of real cases, probably due to the deficiencies and difficulties of testing a significant sample of the whole living population. However, we also remark that COVID-19 mortality data depend on the currently unknown transmission rates, which emphasizes the importance of qualitative agreement to test novel modeling approaches. To this end, we show that it is possible to rescale our simulation results to accurately match the order of magnitude of the reported infected case data Fig. 2 , though we emphasize that this is purely for the purposes of visualization. We further use our model to assess four illustrative reopening scenarios over four months following 27 February 2020: maintenance of restrictions, relaxation of the lockdown everywhere on 3 May 2020 under two different sets of assumptions, and a combination of maintenance of restrictions in Milan and relaxation elsewhere in Lombardy. We still consider the changes in parameter values induced by the sequential restrictions (see Section 3) and the lockdown relaxation is modeled by setting ν s,e,r = 2.175 · 10 −2 km 2 · persons −1 · day −1 and β i,e = 9.0 · 10 −5 persons −1 ·day −1 (scenario A), and β i,e = 6.6 · 10 −5 persons −1 ·day −1 (scenario B). Scenario A is a pessimistic scenario, which assumes that the population contact rate is similar to early-outbreak levels. Scenario B is more optimistic and assumes that the generally greater public awareness of preventative measures (such as mask-wearing and social distancing) translates to greater success in limiting contact, despite increased mobility. Fig. 4 shows the resulting outbreak dynamics for these four reopening scenarios. Our simulations suggest that relaxing the lockdown restrictions in the entire region may cause severe and rapid growth in the Milan area. However, major urban zones far from Milan (e.g., Brescia and Bergamo) just experience a marginal increase in growth and still show a favorable trend in time. Conversely, if we maintain the lockdown restrictions in Milan and relax them elsewhere, the outbreak shows more favorable dynamics, similar to those obtained for Brescia and Bergamo. Thus, our results suggest that maintenance of lockdown measures in high-population, high-density areas like Milan may be necessary for longer times to effectively arrest the spread of contagious diseases like COVID-19. We have introduced a compartmental PDE model describing spatio-temporal propagation of disease contagion and applied it to the 2020 outbreak of COVID-19 in Lombardy specifically. Our simulations are intended to provide a proof-of-concept of the potential of PDEs for regional modeling of the outbreak. Nonetheless, they show good qualitative agreement with reality, accurately predicting the outbreak dynamics in different areas and recreating the transmission path in time and space. We then used the model to examine some possible reopening scenarios, which suggested that reopening may be best determined based on local population and contagion dynamics, and not by a one-size-fits-all approach. Our model is in a very early stage, with ample room for improvement. We plan to consider non-constant model parameters and adaptively update them according to measured data using dataassimilation procedures [4] . Indeed, as more reliable data becomes available, we can further extend parameter calibration to fit data for additional model compartments other than the deceased subgroup. A significant issue when considering complex geographical domains is the quantification of the boundary conditions. Should data not be available, boundary conditions can be defined realistically by including 0D SEIR models describing the fluxes to neighboring regions, similar to what was done for cardiovascular problems (see, e.g., [6, 16, 19] ). Additionally, we used population-dependent diffusion terms, but ideally these could also be affected by geographical features (e.g., rivers, mountains, roadways, and railways) [10] . Non-local effects like the ones modeled by fractional operators can also be included [11] . These considerations may be crucial whenever using the model in larger geographical domains. We would also like to extend our framework into more sophisticated compartmental models including, e.g., hospitalizations, patients in intensive care units, or age and biological sex structures [7, 8] . This would further increase the utility of the model, potentially helping decision-makers to determine the allocation of resources among different areas. The present results clearly pinpoint the current standpoints J o u r n a l P r e -p r o o f Journal Pre-proof of virologists, emphasizing the need of restrictions. Finally, the socio-economical costs of lockdowns are not included here, but could ultimately be incorporated in future quantitative analyses, e.g., aiming at the comprehensive optimization of pandemic-arresting measures. Incubation period of 2019 novel coronavirus (2019-CoV) infections among travelers from Wuhan, China Serial interval of COVID-19 among publicly reported confirmed cases The Ensemble Kalman Filter Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand Multiscale modelling of the circulatory system: a preliminary analysis. Computing and visualization in science Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy Partial differential equations in ecology: spatial interactions and population dynamics Numerical simulation of a suceptible-exposed-infectious spacecontinuous model for the spread of rabies in raccoons across a realistic landscape Modeling the dynamics of novel coronavirus (2019-nCov) with fractional derivative Galerkin methods for a model of population dynamics with nonlinear diffusion A numerical method for spatial diffusion in age-structured populations Postitive RT-PCR test results in patients recovered from COVID-19 Serial interval of novel coronavirus (COVID-19) infections Geometric multiscale modeling of the cardiovascular system, between theory and practice COVID-19 and Italy: what next? The Lancet Partial differential equations in action, from modeling to theory Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries Comparison of the spatiotemporal characteristics of the COVID-19 and SARS outbreaks in mainland China. medRxiv The authors would like to acknowledge the work of Marco Demarziani, Luigi Greco, Isabella Atcha, Kasey Cervantes, Sanne Glastra, Shreya Rana, Stefano Minelli, Chiara Macchello, Simona Petralia, Anita de Franco, and Martina Moschella for their crucial help with data acquisition.