key: cord-0678075-7a8vjga6
authors: Wang, Zhenlin; Teja, Mariana Carrasco; Zhang, Xiaoxuan; Teichert, Gregory; Garikipati, Krishna
title: System inference via field inversion for the spatio-temporal progression of infectious diseases: Studies of COVID-19 in Michigan and Mexico
date: 2021-04-29
journal: nan
DOI: nan
sha: c9273701a25170c5d408a8c7a094ee2d707ca465
doc_id: 678075
cord_uid: 7a8vjga6

We present an approach to studying and predicting the spatio-temporal progression of infectious diseases. We treat the problem by adopting a partial differential equation (PDE) version of the Susceptible, Infected, Recovered, Deceased (SIRD) compartmental model of epidemiology, which is achieved by replacing compartmental populations by their densities. Building on our recent work (Computational Mechanics, 66, 1177, 2020), we replace our earlier use of global polynomial basis functions with those having local support, as epitomized in the finite element method, for the spatial representation of the SIRD parameters. The time dependence is treated by inferring constant parameters over time intervals that coincide with the time step in semi-discrete numerical implementations. In combination, this amounts to a scheme of field inversion of the SIRD parameters over each time step. Applied to data over ten months of 2020 for the pandemic in the US state of Michigan and to all of Mexico, our system inference via field inversion infers spatio-temporally varying PDE SIRD parameters that replicate the progression of the pandemic with high accuracy. It also produces accurate predictions, when compared against data, for a three week period into 2021. Of note is the insight that is suggested on the spatio-temporal variation of infection, recovery and death rates, as well as patterns of the population's mobility revealed by diffusivities of the compartments.

Classical epidemiological models, such as the Susceptible-Infected-Recovered (SIRD) model [1] , are ordinary differential equations (ODEs) defined by specifying the compartmental sub-population numbers over some geographical region. Spatial effects have typically been introduced by resolving smaller regions and treating them individually [2] [3] [4] [5] [6] . During the long-lasting and widespread epidemics, such as the COVID-19 Pandemic, the effects on the infection rate of imposing-and then lifting-mobility restrictions and social distancing mandates revolve on the question of the time and spatially varying mobility of the population. At the finest resolution, this must be approached via agent-based models [7] , using individuals' mobility data. However, this data is not available for the entire population, and contact tracing campaigns face challenges of recruiting workers, access, technology, as well as socio-political resistance. Against these difficulties, an intriguing question to explore is whether simple reaction-diffusion models can detect the evidence of mobility in these data. Such an approach must start with a partial differential equation (PDE) version of the epidemiological models, which is easily defined by converting compartmental sub-populations to densities over sub regions by normalizing with the corresponding areas. To address the mobility of the population, diffusion terms are introduced to the SIRD model, which is transformed to a set of reaction-diffusion PDEs in two spatial dimensions [8] [9] [10] [11] .

The widespread availability of data in the public domain [12] [13] [14] [15] [16] [17] [18] [19] has spurred a widespread interest among computational and data scientists, who have sought to test and refine their methods against these repositories. This has opened up the possibility that advances in computational and data science may contribute to the existing and rapidly expanding body of work in epidemiology, in inferring the dynamics of COVID-19 and making projections. We have similarly sought to build off our recent work in data-driven and machine learning approaches [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] and presented a class of system identification techniques for inference of ODE and PDE forms of the SIRD model, as well as Bayesian neural networks for representation and uncertainty quantification-guided prediction [11] . That work focused on the US state of Michigan. In this communication, we revise our approach for inference of the PDE SIRD model with temporal and spatial evolving parameters and diffusivities. Importantly, instead of global polynomial representations of PDE SIRD parameters over the spatial and temporal domains, we adopt field inversion over time intervals that coincide with the time steps of our underlying numerical implementation. This affords much greater accuracy over the global polynomial ansatz. Adjoint-based gradient optimization for field inversion of parameters at each time step replaces the use of stepwise regression-based system identification in our previous work. We find that the improved accuracy with respect to the data over the time interval of inference, as well as of the predictions, is worth the increased expense. We have brought abundant high-quality, public domain, data [12] [13] [14] [15] [16] [17] [18] [19] on the evolution of COVID-19 in both, the State of Michigan, with a population of 9.98 million, distributed in 83 counties, over 250,493 km 2 , and the country of Mexico, with a population of 126 million, distributed in 32 geographical entities (31 states, plus Mexico City), on 1,972,550 km 2 . The temporal resolution by days and spatial resolution by counties/states have allowed us to study the mobility in these data using our methods of system inference.

In Section 2 we review our previous work of system inference for the spatio-temporal SIRD model first, and then extend it by incorporating temporal and spatial parameters and diffusivities 2 using a finite element representation. The PDE SIRD model-constrained inference are presented in Section 3. Section 4 is on data preparation. The results for inference of classical SIRD parameters as well as the diffusivities, and forward prediction are presented in Section 5. Our conclusions appear in Section 6.

2 Compartmental differential equations models of infectious disease dynamics

We begin with the conventional SIRD compartmental epidemiology model. The population, taken to remain constant at N , is divided into four disjoint compartments with time-dependent subpopulations: S(t) for susceptible, I(t) for infected, R(t) for recovered and D(t) for deceased individuals. The governing system of ordinary differential equations (ODEs) is:

This is the canonical form of the model where the sub-populations are assumed to be well-mixed so that spatial variations can be ignored over the domain of interest. Here β is the infection rate, µ is the recovery rate, γ is the rate of immunity loss, and α is the death rate.

We have extended the SIRD model to a system of partial differential equations (PDEs) in two spatial dimensions using the same compartments [11] . However, the population variables are now replaced with spatio-temporally varying densities, S(x, t), I(x, t), R(x, t), D(x, t) defined as numbers per unit area.

where D S , D I , D R are diffusivities of the corresponding compartments, and represent the mobility of the sub-population via random walks. We define (•) = (•)/ Ω dA where Ω is the domain of study: either the lower peninsula of the State of Michigan, or the territory of the country of Mexico. Furthermore the population constraint holds:

In what follows of this communication, we only consider the PDEs SIRD model, and, for the sake of readability, we dispense with the hats on the compartments.

We adopt the weak form, and specifically, the finite element framework for the above system of PDEs. For a generic, finite-dimensional field u h , the problem is stated as follows: 

where n b is the dimensionality of the function spaces S h and V h , and N a represents the basis functions. To obtain the Galerkin weak forms, we multiply each equation in strong form in (6-9) by a weighting function w h S , w h I , w h R , w h D , respectively, integrate by parts, apply boundary conditions appropriately, and use the Backward Euler method for time-discretization with (•) n denoting a discretized quantity at time t n and ∆t being the time step. See Ref. [11] for details. This leads to:

where the boundary terms vanish because we assume that the sub-populations do not leave the region-zero flux boundary conditions. In our previous work [11] , we have characterized the coefficients to vary via a global-in-time polynomial basis. While the inferred model reproduced the trends, there was a notable error over time of the statewide sub-populations estimates S(t, x), I(t, x), R(t, x), D(t, x) obtained by forward simulation with inferred quantities (See Figure 14 and 15 in [11] ).

Additionally, the highly complex geometry of the State of Michigan, and of Mexico (See maps in Figure 1 ), and potentially highly nonuniform distributions of the coefficients in space makes it challenging to characterize the coefficients with simple basis functions. Global polynomials in space could not sufficiently resolve the emergence and disappearance of "hot spots" and "cold spots" [11] . In this communication, we allow the coefficients β, γ, µ, α, D S , . . . , D R of the PDE SIRD model to vary over space via finite-dimensional, locally supported representations as we do for the primary variables S(t, x), I(t, x), R(t, x), D(t, x). Further more, we allow the coefficients to vary daily, leading to:

where, as for the primary variables, the subscripts (•)n denote the coefficients on day n. With this, the PDE SIRD equations become:

where the PDE SIRD parameters are interpolated from nodal variables as defined in Equations (15) and (16).

The system inference problem is to invert for the quantities β a n , γ a n , µ a n , α a n , D a Sn , . . . , D a Rn . Since these quantities are interpolated via Equations (15) and (16) to be expressed as the corresponding fields β h n , γ h n , µ h n , α h n , D h Sn , . . . , D h Rn in (17) (18) , the system inference problems is one of field inversion. It is stated in Equations (21) (22) as:

Given (15) (16) , at each t n : β a n , . . . , D a Rn np a=1 = arg min (β a n ,...,D a Rn )a=1 np i , such that (17) (18) (19) (20) hold (21) and i is the loss function defined:

where (•) d denotes data for the corresponding quantity. Due to the large differences in the magnitudes of different sub-populations, we choose the weights W S , · · · , W D to be:

.

The weights normalize the sub-populations and prioritize regions with higher infected populations. These regions are of greater interest for studying the progression of the disease as they tend to have a higher population density and, therefore, infected populations. This PDE-constrained optimization problem is solved iteratively, and requires the gradient of the PDE constraint, Equations (17) (18) (19) (20) , with respect to parameters. We adopt classical adjointbased gradient optimization. This approach involves a single linear solution of the adjoint equation of the original PDE constraint at each iteration, followed by solution of the fields to be inverted: (β h n , γ h n , µ h n , α h n , D h Sn , . . . , D h Rn ) and the updated forward solution S h n , I h n , R h n , D h n . In this work we use the L-BFGS-B optimization algorithm from SciPy [30] and the dolfin-adjoint software library [31] to compute the gradient.

First, we constructed two-dimensional meshes for Michigan and Mexico that fully resolve the counties/states as shown in Figure 1 . The data are available as cumulative sub-population numbers I d n , R d n , D d n at the county/state level. We used a uniform density of each sub-population to compute I d n , R d n , D d n within the county/state, and applied Gaussian filtering to smooth the discontinuities at the county/state boundaries. Note that the discrete Gaussian filter can not be applied in a straightforward manner to unstructured meshes. Starting with a field u that represents any of the four sub-population densities, and G(x 0 , x) = 1 2πσ 2 e − ||x|| 2 2σ 2 as the two dimensional Gaussian distribution function centered at 0 with standard deviation σ, which is related to the kernel size in the discrete Gaussian filter, we scale the filtered solution denoted by u(x 0 ) at each finite element node:

The spatio-temporal evolution of these fields was used in the system inference problem as described in Section 3.

7 Figure 2 shows the sub-populations S(x, t), I(x, t), R(x, t), D(x, t) in both Michigan and Mexico obtained by forward simulation with inferred quantities compared with data on December 29, 2020 (t = 281 days), where t = 0 is March 23, 2020, the start of the lockdown in Michigan. The inferred model for Michigan accurately replicates the initial burst of disease and the following multiple waves around Detroit (please see the SI movie: michigan prediction.mp4). It also captures the second burst in the southwest of Michigan around the city of Grand Rapids. The high burden of the disease in these, the largest and second largest cities, respectively, in Michigan, reflects wellknown socio-economic challenges related to Detroit in particular, and more generally reflected in other urban centers. Similarly, Mexico City, with highest population density in Mexico (6, 200/km 2 [18] ), was the worst affected area in that country and dominated the evolution of the disease (See SI movie: mexico prediction.mp4). The low error between the simulation and data leads to greater confidence in the inferred parameters. Figure 3 shows the inferred infection rate, death rate and the recovery rate in Michigan every 70 days starting from t = 0 (the time-resolved dynamics are shown in SI movie: michigan parameter.mp4). The evolution of these inferred parameters reveals that the population's infection rate, β(t), declined from the initially higher values in highly infected areas (such as Detroit), and spread to the western parts of Michigan. The death rate was mostly stable after May 2020 (t > 69), and remained low in the more highly infected areas. This can be attributed to the ramp up of the public health campaign, hospitalizations and emergency response of the medical system, and prioritization to the more highly infected areas. The recovery rate around Detroit city evolved in multiple stages: increasing→ decreasing → increasing, which was consistent with the multiple waves reflected in the data on the recovered population in this region (SI movie: michigan prediction.mp4).

At the finest resolution, the mobility of the population during disease evolution may be approached via agent-based models refined to resolve individuals. However, given the difficulties encountered in effective contact tracing, and its acceptance by the population [9, 32, 33] , an intriguing question to explore is whether simple reaction-diffusion models can detect the evidence of mobility in these data. Figure 4 shows the inferred diffusivities of the susceptible, infected, and recovered sub-populations. Note that for field inversion, the population density data for each compartment was taken to be uniform within each county/state, since no finer grained information was available, and then subject to Gaussian smoothing before inference. Thus the density gradients, which drive the inference of diffusivities, arise at the counties/states scale more than they do at the intra-county/intra-state. Accordingly, the inferred diffusivities are meaningful on this scale. The lower Peninsula of Michigan is about 446 km long from north to south and 314 km wide from east to west-scales that can help place the diffusivities in Figure 4 in perspective. The mobility of the infected population was always high around the highly infected areas. In Michigan, this infected population gradually shifted to the southwestern part of the state from the initial burst around Detroit. This finding is consistent with the second burst around Grands during the evolution of the pandemic. The recovered population demonstrated a similar pattern of mobility, and was more active in the southern part of Michigan around the more highly infected regions. On the other hand, susceptible population closely tracks the total population. Since the population at large has low mobility, the susceptible population's mobility is low in high population density areas. See SI movie michigan prediction.mp4 for these dynamics. Figure 5 show the inferred infection rate, death rate and the recovery rate of the inferred model for Mexico. We can clearly see the spreading of the disease from Mexico City. Similar to the case of Michigan, infection rates, and to a lesser extent, death rates, were relatively lower in Mexico City, which is the most densely populated region of the country, than that in the surrounding cities. The recovery rate was high in Mexico city due to the relatively greater resources of the medical system there. The infection and death rates tended to be stable for five months following March 23, 2020, and the recovery rate gradually increased in more areas. Notably, far from the Mexico city, Baja California also displayed a high inferred rate of infection. We suspect this to be because it borders California, USA, and the international border restrictions did not contain the spread of the virus between the two regions. Unlike Mexico City, the death rate remained high, and the recovery rate did not increase to levels comparable to the capital, perhaps because of the looser restrictions in this popular tourist destination.

The diffusivities of the corresponding sub-populations of the inferred model for Mexico are shown in Figure 6 . Similar to Michigan, the mobility of infected and recovered sub-populations are higher around the highly infected Mexico City. Mexico is about 3000 km long from north to south and 1900 km wide from east to west-scales that can help place the diffusivities in Figure 6 in perspective. Unlike the case in Michigan where there were multiple bursts in different cities, the mobilities of all sub-populations became stable after about 5 months from March 23, 2020. This may reflect differences in the proclivity toward domestic/local mobility of the populations of Michigan and Mexico-two regions with strongly contrasting social, economic and cultural characteristics.

Finally, taking the inferred parameters on the last day used for inference (Day 281), we predicted the evolution of sub-populations for three weeks (Days 282 to 303) using the inferred model. Figures  7 and 8 show the predicted spatio-temporal evolution of the infected population against the raw data for both Michigan and Mexico. The inferred models captured closely the evolution of the infectedsub-populations, indicating that the dynamics of the disease tended to be steady in January 2021. The prediction of recovered and deceased sub-populations are shown in Figures 9 and 10 under Supplementary Information. 

This communication builds upon our previous work [11] on system inference and machine learning from data to study the progression of COVID-19 across the state of Michigan. We extended the PDE SIRD model by allowing the infection rate, death rate and the recovery rate, as well as the diffusivities of the susceptible, infected, and recovered sub-populations to vary over space and time. Using field inversion to infer the parameters as finite-dimensional fields on time scales of a single day, we obtained models to predict the evolution of disease with high accuracy. This provides us with the ability to analyze the dynamics of the disease through the inferred parameters, and make accurate predictions within a reasonable time frame. Particularly, we can detect the evidence of time and spatially varying mobility of the population through the simple diffusion-reaction models instead of the relying on the agent-based models which require individual's mobility data. The latter can prove challenging, technically as well as politically, to obtain.

As discussed in Section 5, our inferred models capture the geographical spread of infection, the number of deaths and the size of the recovered population starting from one highly infected area to its surrounding cities and eventually spreading to further areas. Particularly, the higher infection and death rates in areas with low infection at later times suggests that more attention is needed in such locations. This may be due to a lack of medical services, or a lack of compliance with mitigation strategies. Our inferred models also reveals higher mobility surrounding the highly infected areas suggesting the importance of quarantine and social distancing.

Finite-dimensional representation allows the parameters to accurately capture the spatial dependence, however the non-parametric representation makes the projection of these parameters beyond the data range extremely challenging. A prediction cannot be made with confidence if the dynamics of the disease reflected by these parameters are not stable. Of course, extrapolation is challenging in almost all data-driven methods. One possible alternate is to develope surrogate models of these parameters via time dependent neural networks under the constraints of the SIRD model to learn the spatial variation in time, and thus to make reasonable prediction of the dynamics in the evolution the disease, such as we have demonstrated previously [11] . Nevertheless, without including factors such as mobility restrictions or other mandates, only short time predictions may be accurate. 

A contribution to the mathematical theory of epidemics

Forecasting and uncertainty in modeling the 2014-2015 ebola epidemic in west africa

Examining rainfall and cholera dynamics in haiti using statistical anddynamic modeling approaches

Quantifying the impact of human mobility on malaria

Modeling the worldwide spread of pandemic influenza: Baseline case and containment interventions

The mathematics of infectious diseases

A taxonomy for agent-based models in human infectious disease epidemiology

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to covid-19, mathematical analysis, and numerical study

An agent-based computational framework for simulation of global pandemic and social response on planet x

Cross-diffusion-induced patterns in an sir epidemic model on complex networks

System inference for the spatio-temporal evolution of infectious diseases: Michigan in the time of COVID-19

Covidnet: To bring data transparency in the era of covid-19

Michigan State Coronavirus Data

Latest Map and Case Count -The New York Times

The Institute for Health Metrics and Evaluation. COVID-19 Projections

Censo de población y vivienda

Conacyt: Covid-19 méxico

Variational system identification of the partial differential equations governing the physics of pattern-formation: Inference under varying fidelity and noise

A perspective on regression and bayesian approaches for system identification of pattern formation dynamics

Variational system identification of the partial differential equations governing microstructure evolution in materials: Inference over sparse and spatially unrelated data

An inverse modelling study on the local volume changes during early morphoelastic growth of the fetal human brain

Discovery of deformation mechanisms and constitutive response of soft material surrogates of biological tissue by full-field characterization and data-driven variational system identification

Machine learning materials physics: Surrogate optimization and multi-fidelity algorithms predict precipitate morphology in an alternative to phase field dynamics

Machine learning materials physics: Integrable deep neural networks enable scale bridging by learning free energy functions

Scale bridging materials physics: Active learning workflows and integrable deep neural networks for free energy function representations in alloys

Machine learning materials physics: Multi-resolution neural networks learn the free energy and nonlinear elastic response of evolving microstructures

Bayesian neural networks for weak solution of pdes with uncertainty quantification

Paul van Mulbregt, and SciPy 1. 0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python

dolfin-adjoint 2018.1: automated adjoints for fenics and firedrake

Privacy concerns can explain unwillingness to download and use contact tracing apps when covid-19 concerns are high

Impact of delays on effectiveness of contact tracing strategies for covid-19: a modelling study

We acknowledge the support of Defense Advanced Research Projects Agency (DARPA) under Agreement No. HR0011199002, "Artificial Intelligence guided multi-scale multi-physics framework for discovering complex emergent materials phenomena"