Accurate modeling of spatiotemporal spread is critical when faced with an epidemic at a global scale such as the recent outbreak of COVID-19. The disease is caused by the SARS-CoV-2 virus, whose outbreak was first witnessed in the city of Wuhan in the Hubei province of China in December 2019. Since then the virus has spread rapidly to other parts of the world over a span of just a few months, thereby causing all countries to respond by imposing strict gathering restrictions, travel bans, stay-at-home orders and occasionally even curfews to counteract the virus spread. Given the fast spread and limited data collection, it has become necessary to advance our responding capabilities in various fields like contact tracing, medicine, healthcare, genome tracking and spread modeling.
Epidemic spread has been traditionally modeled via compartmentalized models, e.g., the SIR model [14] . However these models do not account for under-reporting of confirmed cases or deviations from natural evolution of the disease caused by various intervention policies imposed by humans to counteract the spread. Estimating the model parameters is generally done via statistical analysis which assumes free evolution of the disease in the absence of any human intervention and requires guessing the initial hidden states of all compartments [13] . In this work, we present our PolSIRD model which leverages intervention and disease spread data to model the disease spread. We account for under-reporting of cases via an observation mechanism and allow the model to learn the initial hidden states of all compartments along with other parameters via gradient-based training. Lastly, PolSIRD accepts intervention policy inputs and learns their effects on disease spread parameters. We apply our model to the spread of COVID-19 and also provide actionable guidance on intervention policies via counterfactual simulations from our model.
Modeling the spatiotemporal spread of diseases has been studied in epidemiology with the longest standing models being variants of the compartmentalized SIR model proposed by Kermack and McKendrick [14] . The SIR model is a continuous-time markov model useful for studying the evolution of a disease in a single homogeneously mixing population. Over time, variants of the SIR model have been proposed and analyzed to account for discrete time, latent stages of the disease (SEIR) and lack of immunity after recovering (SIS) [1] . More advances resulting from modeling the global spread of Influenza have additionally accounted for the spread of disease via transportation between multiple heterogeneously mixing populations [23, 20, 12] . While the parameters of these models were acquired by statistical analysis of the disease spread, we learn all parameters of PolSIRD end-to-end via gradient-based training on observed data.
Other approaches use stochastic modeling approaches based on point processes (e.g. Hawkes processes) to model disease spread in multiple populations spread across geographic locations [16] , however these require availability of event-level data which is often not be available due to underreporting of cases. Other point process based approaches use aggregated event data but do not model all compartments, nor the effects of finite population sizes [19] . While modifications to Hawkes processes have been proposed to have a Hawkes-SIR model [22] , these only roughly approximate the desired compartments over a finite population and have essentially been used to estimate the diameter of disease spread on the underlying graph of geographic locations. The more recently proposed GLEaM model uses a Multinomial distribution to exactly model a finite population, however, unlike PolSIRD, GLEaM does not account for under-reporting and uses statistical estimates of spread parameters instead of learning them end-to-end from data. We refer the interested reader to recent epidemiological reviews [24, 8, 27] which survey the advances in this field in more detail.
More recently, the outbreak of COVID-19 has caused numerous works to appear focusing on statistically estimating the spread and death rates from Wuhan data [26] and on using the SEIR model (SIR augmented with latent stage of disease) to estimate the spread in Wuhan and the rest of China [28, 6] . Other works have used point process based models to estimate the spread locally [21] . However, the latent stage of COVID-19 is not easily observable since a latent individual is asymptomatic and yet infectious within local populations or between populations via travel. Hence, COVID-19 is plagued by under-reporting of confirmed cases which can be as high as 86%, as was the case in Wuhan [18] . This along with the non-availability of enough testing equipment has made COVID-19 spread across the whole world in a matter of three months. While the recent SuEIR model [29] has modified the SEIR model to account for under-reporting of cases, it cannot simultaneously learn about intervention policies.
Many works focus exclusively on studying effects of a particular intervention policy like restrictions on human mobility and foreign travel [5, 7] and quarantine [9] in various countries. While the recent DELPHI model [17] also accounts for an intervention policy, it does not provide any model for combining multiple intervention policies. Our work differs from these works in that it assumes the presence of multiple intervention policies and expects only a fraction of cases to be reported. PolSIRD learns the spread, policy and reporting parameters end-to-end directly from observed data. Hence, after learning our model we can estimate the contribution of individual intervention policies towards curbing the COVID-19 spread in the United States. We are also able to run counterfactual simulations and make predictions on the effects of lifting certain intervention policies in the near future as is planned for many US states. The SIR compartmental model [14] models the evolution of a disease in a fixed-size population in continuous time using ordinary differential equations. We describe here a discrete-time variant proposed later [20, 1] . Assuming a population size N , the SIR model maintains three population compartments at all times: (a) S: Number of people susceptible to the disease, (b) I: Number of people infected by the disease and actively spreading it further, and (c) R: Number of people who have been removed either due to death or due to a full recovery leading to immunization. The state transition for individuals follows from S → I → R, hence the name SIR. At any time t, an individual can be in exactly one compartment and the total population is always conserved,
The Markov evolution equations of the SIR model with time step ∆t are:
where the first equation assumes uniform infectious contact between individuals in compartments S and I and the last equation assumes uniform removal (via recovery or death). The parameter λ ≥ 0 governs the infection rate and the parameter γ ∈ [0, 1] is the removal rate. Since epidemic reported cases are generally recorded on a day-to-day basis, we will use ∆t = 1 without loss of generality.
Under the SIR model, the average length of a person's infection is given by 1/γ. Another important quantity is the basic reproduction number of the disease given by r 0 = λ/γ, i.e. the average number of people a single infectious person will infect over the course of their infection. The number r 0 decides if the disease will persist in the population or die out. In general, r 0 < 1 implies the disease automatically getting extinct, r 0 = 1 makes the disease endemic thereby stably persisting in the population and r 0 > 1 makes the disease grow exponentially and become an epidemic. 
S Ĩ I R D C P ol {λ,λ} β I γ γ γ(1 − γ C )β D w, b
Epidemic spread data is generally plagued by under-reporting due to lack of reliable testing technology and/or due to availability of sufficient testing kits. While traditional SIR models do not account for under-reporting, we assume that the confirmed counts are only a fraction of the true total cases. We keep the susceptible (S), unreported but actively infectious (I) and the removed (R) compartments. Apart from these, we add two new compartments to account, namely the reported infectious casesĨ and the reported deathsD. However, while infected individuals can recover or die and leave theĨ compartment, the actual reported numbers do not track the active infectious people but rather the total number of individuals ever confirmed to be infectious till date. Hence, we also add an observableC to keep track of the total number of confirmed cases. Note that we callC an observable and not a compartment since individuals in it belong toĨ,D and R compartments and it double counts them for observation purposes.
Temporal evolution: The population conservation can now be written as:
Note that the confirmed cases observable (C) is not part of the population conservation equation since it logs individuals already present in other compartments. The Markov evolution equations of the system with state (S, I,Ĩ, R,D,C) for time step = 1 day are given as:
) which are also described graphically in figure 1. As eq 6 shows, we now model infections to susceptible people from both I andĨ compartments with their individual infection rates λ,λ ≥ 0 respectively. While susceptible people go first into the unreported infectious compartment I, they can eventually get tested and go to the reported infectious compartmentĨ with a testing fraction β I ∈ [0, 1] (eqs 7, 8) . Note that the I compartment also contains people who have been exposed and are infectious but not yet symptomatic, e.g., as found for the latent stage of COVID-19. Individuals in the infected compartments I andĨ are removed (by recovery or death) into the R compartment with a rate γ ∈ [0, 1] (eqs 9). However, a fraction of the deaths from confirmed infected casesĨ are also reported in theD compartment. Assuming a cure rate of γ C ∈ [0, 1] and a death reporting rate of β D ∈ [0, 1], the number of reported deaths per day are given as γ(1 − γ C )β DĨ (t) (eq 10). Note that while we have both γ C and β D as learnable parameters, one needs to observe infections, recoveries and deaths to be able to learn the parameters separately. However, since only the confirmed infections C and confirmed deathsD are being reported reliably for COVID-19, one cannot disambiguate the effects of γ C and β D and our model learns the product (1 − γ C )β D jointly as a parameter.
Influence of intervention policies: While most existing models simulate free evolution of a disease, the data collected generally comes from an initial freely evolving disease followed by a subsequent period of controlled spread due to execution of intervention policies by humans. This requires us to model the influence of intervention policies, e.g. executing stay-at-home orders or closing gyms and movie theaters. Since the most rapid spread is caused by community spread, most policies are targeted towards curbing the infection rates λ andλ within a population. Denoting the set of all policies under consideration by P, we represent the time at which policy p ∈ P was enacted as τ p . A single policy p is represented by four learnable parameters: w p ,w p , b p ,b p where w p (w p ) represents the decay rate and b p (b p ) represents the steady-state reduction for λ(λ) due to the policy p. While different policies can interact with each other in complex ways, we approximate the effect of multiple policies with a first-order approximation which assumes that all policies independently reduce the spread rates and multiplies the reduction to λ (orλ) from all the enacted policies. Given that all policies under consideration are generally enacted within a short time frame, as has been the case with COVID-19 in the US, it is hard to disentangle the individual effects of policies under any model which considers combining multiple policies. This becomes even harder if more complex forms of interactions are modeled between policies, hence we leave the modeling of more complex policy interactions for future work and stick to our simple and more interpretable model in this work. With multiple policies being enacted at a node, both λ andλ become decaying functions of time, given by:
where sgn(p) is +1 if the policy p is being enacted and −1 if a pre-enacted policy p is being lifted.
Note that the above form ensures that enacting any policy always reduces the value of λ andλ vice-versa when lifting policies, but more importantly in a specific way. We describe the effect of eq 12 for a single policy on λ and the generalization with multiple policies automatically follows:
• On being enacted, an intervention policy generally causes some immediate spread reduction due to a fraction of the population adopting the policy immediately. In our model, the coefficient of λ at time t = −∞ is 1, thereby causing no reduction to the natural rate. Since the decision to enact a policy is generally rapidly taken in a few days, τ p is not known much ahead of time before the policy is about to be enacted. Hence, we keep t at −∞ till τ p . When the policy is enacted at time τ p , plugging in t = τ p immediately drops the coefficient of λ to 1+bp 2 . • However, a policy continues to further reduce the spread slowly over time as the remaining population adopts the policy. Our model causes such a slow decay after time τ p and asymptotically decays the coefficient of λ to its steady state value of b p with an exponential rate given by w p .
• Lastly, one can perform counterfactual simulation by lifting policy p after λ has settled to its steady state value after p being enacted. This follows a similar trend of a certain fraction of people adopting the lifting immediately while other more cautious people slowly adopt it after observing the effects of lifting. On lifting policy p, our model increases the coefficient of λ suddenly from b p to 1+bp 2 and then asymptotically increases it back to 1 as time passes.
We use λ(t) andλ(t) from eqs 12,13 in the temporal evolution eqs 6-11 while unrolling over time.
Initialization of PolSIRD: Before unrolling the recurrent PolSIRD model over time, one needs to guess an initial estimate for the number of people in each compartment. This is a particularly challenging task since the true total number of people in any compartment are never observed at any time. While existing approaches [3, 29] have chosen suitable initializations via approximate searches and hyperparameter tuning, we take a different route and take the initial values in each compartment at t = 0 as learnable parameters, i.e., our model learns the fractions s 0 , i 0 ,ĩ 0 , r 0 andd 0 such that they sum upto 1. While unrolling the model from t = 0, we set the initial state of the model with total population N as:
where the initialization forC(0) is given by adding all reported active infections (Nĩ 0 ) to all the reported cases which left the compartmentĨ to go toD or to R. Since we know that people leave fromĨ at a rate γ, a fraction of which entersD with rate γ(1 − γ C )β D , having the total reported deaths at time t = 0 as Nd 0 implies that the total number of people who ever leftĨ =
(1−γ C )β D .
The PolSIRD model takes as input the time-steps at which each policy was enacted and unrolls the temporal evolution model defined by eqns 6-13 for T time steps. The model is trained to learn all its parameters θ = {λ,λ, γ, β I , β D , γ C } ∪ {w p , b p ,w p ,b p } p∈P ∪ {s 0 , i 0 ,ĩ 0 , r 0 ,d 0 } by minimizing the Mean Absolute Error Ratio (MAER) on the number of confirmed casesC(t) and the number of reported deathsD(t) at all time steps:
where the max(·) term in the denominator serves to handle nodes with 0 reported cases and a hat(ˆ) on the top denotes the observed ground truth value. While one can train the model independently at any location, in practice we train models jointly across multiple locations (e.g. on all US states). While each location maintains its own copy of the temporal evolution parameters and the initial state parameters for compartments since these can vary significantly across locations, we enforce that policies have similar effects at all locations and leverage additional data for learning about policy effects by sharing the same copy of policy decay rates and steady-state reductions amongst all locations for each policy 3 .
We apply our PolSIRD model to predicting the spread of COVID-19 and the corresponding spread parameters for all states in the United States. We use the following datasets for training our model: Comparing with baselines: We compare our model specifically against several baselines generated by other institutes and being used to inform the decisions made by the Centers for Disease Control and Prevention (CDC): (a) SuEIR model proposed by Zou et al. [29] , (b) GLEaM model proposed by Balcan et al. [2] and, (c) DELPHI model proposed by Li [17] . Since the CDC only maintains predictions of future deaths, we directly use the predicted death cases from all the baselines being used by the CDC (taken from [11] ) and compare them to our model's predictions. For this we only use the predictions made from May 2nd onwards and ensure that we use versions of all models which have only been trained on the time-series uptil May 1st. We show the MAER and the Root Mean Square Error (RMSE) for all models in Table 1 averaged across all the US states. We also additionally show both metrics averaged across the best 40 states for each model. This removes certain states which have abnormalties in reporting or very few reported deaths causing the models to have a high error on those states. Our PolSIRD model clearly outperforms all the other state-of-the-art baselines in terms of both metrics on all 51 states (counting the District of Columbia as a state) and also on the top-40 states. We show the death predictions from all models on several states in figure 2. 
We report the parameter values learnt by our model averaged over the US states in Table 2 . Our estimated mean disease lifetime given by 1/γ = 17.86 days is a lies in the statistical range of 17.8 to 24.7 days estimated by [26] . Based on β I , our model estimates about 61% reported confirmed [13] indicating a higher spread rate in the absence of any intervention. Since both β I and r 0 values vary considerably across states, we visualize them across all US states in figures 3a and 3b respectively. The high values of r 0 in New York, New Jersey, Texas, Atlanta and Wisconsin are consistent with the known high spread rates in these states. While the reporting fractions are high in Washington, Louisiana, Illinois and New York, many other states suffer from heavy under-reporting. A first actionable measure for such states would be to provision more resources for testing and reporting cases in order to gauge the spread of COVID-19 more accurately. 
Effects of policies: We next analyze the effects of intervention policies on the infection rates and the reproduction number. Table 3 summarizes the steady-state decays to λ,λ from each policy. We observed from our experiments disambiguating the contribution of various policies is hard since they were enacted between 6-20 March with very short durations between them to observe their individual effects. Hence, we focus on the overall steady-state decay due to all policies combined together in Table 4 . Overall, the policies have been quite effective leading to a reduction in spread of about 79% from reported cases and 93% from unreported cases. This is expected because people with confirmed infection tend to isolate themselves while asymptomatic unreported cases unknowingly spread the disease. This unreported spread gets curbed to a larger extent by intervention policies. The new reproduction number r 0 has dropped down to 1.159 averaged across the US states, which while much smaller than its natural value is still > 1, indicating that the disease continues to progress, albeit at a slower rate. However several states now have an r 0 value less than 1 (see figure 3c ) and this has led many state governments to begin formulating plans to lift the policies and re-open the states.
Counterfactual re-opening of states: Most re-opening plans are multi-stage and remove intervention policies sequentially. At the time of writing this manuscript many states are already in stages 1 or 2 of re-opening which prep workplaces to follow social distancing guidelines and allow opening small low-risk workplaces and businesses respectively [25] . To understand the effect of eventually moving into stages 3 (opening restaurants for dine-in, gyms and entertainment venues) and 4 (lifting stay-athome order, bans on large gatherings and public school closures), we now perform counterfactual simulations from our trained model to lift the intervention policies. We execute three counterfactual stage 3 and 4 plans as follows (all stage plans are in reference to May 1, 2020): Plan 1: Enter stage 3 in 60 days and stage 4 in 120 days, Plan 2: Enter stage 3 in 90 days and stage 4 in 120 days, and Plan 3: Enter stage 3 in 90 days and stage 4 in 150 days. The results of our counterfactual simulations when all policies are in place and when the above plans are followed are shown on several states in figure 4 . While the disease progression is slowed under policy intervention, we observe that the number of infected cases begins to rise exponentially if the policies are removed following any of the above three plans; the key takeaway being that it could be detrimental to remove the intervention policies any time in the near future. To be cautious and not cause a second wave of COVID-19 infections, we recommend deferring the stages 3 and 4 of re-opening to more than 150 days (from May 1, 2020) in at least the states which have a current r 0 ≥ 1 (figure 3c). 
We have presented a mathematical model, PolSIRD, which models the spread of epidemics within a population. Our approach accounts for reporting mechanisms, consequent under-reporting and the effects of enacting intervention policies. We successfully outperform the methods currently used by the Center for Disease Control and Prevention for making spread predictions on COVID-19 and provide actionable guidance based on counterfactual simulations on lifting of intervention policies.
Our work broadly aims at studying the spread of COVID-19 in the United States. We model the spread of the disease under intervention policies and under-reporting while trying to provide actionable guidance on the recent plans being formulated to re-open the US states in the next 3-4 months. We believe that both the state governments and the decision making authorities (e.g., Centers for Disease Control and Prevention) will potentially benefit from our work. While we cannot guarantee that all the parameters we have estimated from our model can be deemed exactly accurate, we have made a significant effort to verify that they are in agreement with parameters being estimated from other statistical studies. However, they may still be affected by: (a) biases in data reporting procedures, (b) large variations in availability of testing equipment across the US states, and (c) the short time-window within which all policies were enacted in all the US states.