The classical infectious disease model, the SEIR model (Hethcote 2000) , is a variation on the basic SIR model (Anderson and May 1992) . It assumes that all individuals in the population can be categorized into one of the four compartments: Susceptible, Exposed, Infected and Removed, during the period of pandemic. The model describes the evolution of the compartmental populations in time by a system of nonlinear ordinary differential equations (ODE):
The total population S + E + I + R is invariant in time, which we shall normalize to 1 or 100 % in the rest of this paper. Clearly, SEIR is a simplistic temporal model of a given region or country. However, the infectious disease data often provides not only temporal but also spatial information as in the case of COVID-19, see (Dong, Du, and Gardner 2020) . A natural idea is to elevate SEIR model to a spatio-temporal model so that it can be trained from the currently reported data and make more accurate real-time prediction. See (Roosa et al. 2020) for temporal modeling on cumulative cases of China and real-time prediction.
In this paper, we set out to model the latent effect of inflow cases from the geographical neighbors to capture spacial spreading effect of infectious disease. For the practical reason that the inflow data is not observable, machine learning methods such as regression and neural network are more suitable. As widely adopted in time-series prediction problem, linear statistical models, such auto-regressive model (AR) and its variants are standard methods to forecast timeseries data with some distribution assumptions on the time series. And the Long Short Term Memory neural networks model (LSTM) (Hochreiter and Schmidhuber 1997) for the natural language processing problem, can be applied to time series data, especially disease data. With additional spatial information, the graph-structured LSTM models show a better performance on spatio-temporal data. See the application to influenza data (Li et al. 2019 ) (Deng et al. 2019 ) and crime and traffic data (Wang et al. , 2018 . However, such neural network models have a demand for a large training data to optimize the high dimensional parameters. Yet the reliable daily data of COVID-19 in the US begins after March 2020 and limits the temporal resolution. Applying space-time LSTM models (Li et al. 2019; Wang et al. 2018) directly to COVID-19 may lead to overfitting. In light of the shortage of data of COVID-19, we shall derive a hybrid SEIR-LSTM model with much fewer parameters than space-time LSTMs (Lai et al. 2017 ) (Wu et al. 2018 ).
In (Yang, Santillana, and Kou 2015) , ARGO (AutoRegression with Google search trends), a variant of AR, uses the google search trends to generate external feature of ARGO and forecasts influenza data from Centers for Disease Con-arXiv:2010.09077v1 [cs.
LG] 18 Oct 2020 trol of U.S.(CDC). ARGO is a linear statistical model that combines historical observations and external features. The prediction of influenza activity level is given by:
whereŷ t is the predicted value at time t, and the optimization part of ARGO is:
where α = (α 1 , · · · , α 52 ) and β = (β 1 , · · · , β 100 ). y t−j (1 ≤ j ≤ 52) are historical values of past 52 weeks and X i,t (1 ≤ i ≤ 100) are the google search trend features at time t. The feature are generated by top 100 of most related trends to influenza from google search at each time. The additional regularization terms to linear regression model helps ARGO optimize. The numerical experiment from (Yang, Santillana, and Kou 2015) shows a better performance than machine learning models such as LSTM, AR, and ARIMA.
The (Li et al. 2019 ) introduces a graph structured recurrent neural network (GSRNN) to further improve the forecasting accuracy of CDC influenza activity level data. From CDC data, the USA is divided into 10 Health and Human Services (HHS) regions to report influenza activity level. These 10 regions are described as a graph in GSRNN with nodes v 1 , · · · , v 10 and a collection of edges based on geographic neighbor relationship (i.e. E = {(v i , v j )|v i , v j are adjacent}, E is the set of all edges). By comparing the average record of activity levels, the 10 HHS region nodes are divided into two groups by relatively active level, the high active group H, and low inactive group L. The two group leads to 3 types of edges between them, L − L, H − L, and H − H, where each edge type has a customized RNN, called edge-RNN, to generate the edge features. There are also two kinds of RNNs for each node group to combine the edge feature with historical values and output the final prediction. Suppose a fixed node v ∈ H. The edge feature of v at time t are e t v,H and e t v,L , which are generated by the average of historical values of neighbor nodes of v in corresponding groups. The edge features are the input of the corresponding edge-RNN of each edge:
Then, the outputs of edge-RNNs are fed into the node-RNN of group H together with the node feature of v at time t, denoted as v t , to output the prediction of the activity level of node v at time t + 1, or y t+1 v :
Our approach: IeRNN model
We propose a novel hybrid spatio-temporal model, named IeRNN, by combining LSTM (Hochreiter and Schmidhuber 1997) and I-equation on a graph structure. The I-equation is a discrete in time model derived from SEIR differential equations. It resembles a nonlinear regression model of time series. The LSTM framework is applied to model the latent geographical inflow of infections. Our IeRNN model, comparing to (Li et al. 2019; Wang et al. 2019 Wang et al. , 2018 , is much more compact.
As a variation to SEIR model, we shall construct additional features I e and E e that reveal the inflow population of infectious and exposed individuals from neighboring regions. Then we augment the SEIR differential equations with I e and E e as:
It still follows that
by normalizing compartmentalized populations to percentages of total population. From (1) and (4), we have
Substituting (6), (7) and (5) in (3), we have a closed Iequation:
The above derivation holds for time dependent coefficients
Let E e = τ I e , and write σ 2 τ as σ 2 . By the explicit Euler and (P + 1)-term Riemann sum approximation, we have a discrete time recursion:
gives the I-model:
If I e ≡ 0 in I-model (10) Generate edge feature with I e The spatial information based on US states map (Andrew 2005) , see Fig. 2 , is formulated as an adjacent matrix G = (g i,j ). If two states v i , v j are neighbors to each other, then g i,j = 1 otherwise is zero. With the variables of graph information, we can define the edge feature of state v i at time t:
where I j,t is the infectious population percentage in state v j at time t.
Then we design an edge-RNN composed of stacked LSTM cells Fig. 3 with a following dense layer Fig. 4 to output I e . The edge feature f i,t is the input of the edge-RNN. The integrated procedure to generate edge feature is illustrated by Fig. 5 , taking California as example. Our model, 
During an epidemic, the rate of infection could change as governments start responding to the epidemic. The infectious rate would start decreasing due to the restrictive policy (partial or full lock-down) being put in place. We model the policy response by changing the parameter, β 1 , in the ODE set (10). By multiplying a control factor β decay (called test decay) to β 1 , we can control the infection levels in the future resulting from various degrees of opening policies.
The policy response β 1 is a function of time to reflect no measure, restricting mass gatherings, reopening, lock-down for different times:
where parameters a, b, c, t 0 , s are learned from fitting historical data. 
We use the COVID-19 data in United States (Dong, Du, and Gardner 2020) to evaluate our IeRNN model in training and testing. From (Dong, Du, and Gardner 2020), we find the state level infectious data in the US. Due to the incomplete recovered cases of US, we use the difference of cumulative cases in each state as the daily infected population. Then we use the population of US (World-Population-Review 2020) to calculate the infectious rate of each state, where we assume the population of a state is constant during the period we concerned with. The data is split into training set (133 days) and testing set (35 days) for model evaluation.
The loss function for training is mean squared error (MSE) of the output of model and true data value:
where the output of modelŷ t has the form (adapted from (10)):
where we have parameters (α, β 1 , β 2 , γ, σ 1 , σ 2 ). Due to the interpretations of SEIR model, these parameter values should range in the interval [0, 1]. We use gradient descent optimizer, Adam (Kingma and Ba 2015) , to train our IeRNN model. In each step, we update the weight of neural networks model (10) and the parameters of loss function (13) separately with different length of step and regularization norms.
To assess the performance of our model, we design a series of numerical experiments to compare the IeRNN with I-equation, temporal LSTM and ARIMA.
Regarding model size, the IeRNN and LSTM have about 4240 parameters while the I-equation and ARIMA have 5 parameters.
Model robustness in training is an important attribute, so that the model performance is not sensitive to initialization of parameters (α,β 1 ,β 2 ,γ,σ 1 ,σ 2 ) during training. We find that the I-equation (I-model with I e = 0) is not easy to learn in the sense that a sub-optimal local minimum is often reached by gradient descent during optimization. With coupling to RNN (I e = 0) in IeRNN, the landscape of loss function is regularized so that a local minimum from any random initialization gives a robust and accurate fit. Fig. 8 shows that I-equation is much less accurate in 1-day ahead prediction than IeRNN. Fig. 9 illustrates the same outcome in 7-day ahead prediction.
In further experiment, we train and test IeRNN and I-equation with randomly initialized parameters (α,β 1 ,β 2 ,γ,σ 1 ,σ 2 ). By repeating the training and testing procedure for 20 times, we compare the average MSE loss for both models. The results in Tables 1 and 2 show that IeRNN performs better for both training loss and testing loss in 1day ahead and 7-day ahead predictions.
I-equation
California training 7.63e-09 8.49e-08 testing 1.26e-08 9.68e-07
Florida training 4.24e-08 3.45e-06 testing 3.59e-08 3.97e-05
Virginia training 3.70e-09 2.60e-08 testing 6.90e-09 1.56e-07 
We compare IeRNN (with β 1 (t)), IeRNN, LSTM and ARIMA on 1-day ahead prediction. IeRNN achieves lower MSE error than LSTM and ARIMA on test set. With policy response function β 1 (t), IeRNN gives further improvement beyond IeRNN with constant β 1 , see Table 3 .
IeRNN LSTM ARIMA β 1 (t) California 1.83e-09 2.45e-09 5.00e-09 1.44e-08 Florida 6.13e-09 7.55e-09 4.68e-08 4.11e-08 Virginia 1.27e-09 1.29e-09 3.37e-09 3.74e-09 
Motivated by weekly forecasting from CDC, we study the 7-day ahead prediction task. The loss function is modified by replacing the I-model by a 7-day delayed version below: Figure 6 : Effect of test decay (policy response multiplier) in test period of 1-day ahead prediction task. The IeRNN is trained through March 3, 2020 to July 14, 2020. The vertical axis is fraction of newly infected people in the population. The horizontal axis is time in unit of days. Figure 7 : Effect of test decay (policy response multiplier) in test period of 7-day ahead prediction task. The IeRNN is trained through March 3, 2020 to July 14, 2020. The vertical axis is fraction of newly infected people in the population. The horizontal axis is time in unit of days.Î
where the output value at time t is influenced by the feature vector I e at time t − 7 and later. With a similar modification of loss function, we adapt LSTM to the 7-day ahead prediction. Table 4 compares IeRNN and LSTM in terms of MSE on testing data.
To study the effect of policy response in IeRNN model on testing data, we multiply the learned β 1 (t) by a constant factor (called test decay) during testing. Fig. 6 and Fig. 7 show the impact to model prediction on test data by adjusting test decay which could control the future trend of infection.
IeRNN LSTM β 1 (t) California 6.79e-09 9.84e-09 1.49e-08 Florida 4.34e-08 4.47e-08 5.74e-08 Virginia 1.16e-09 1.55e-09 1.54e-08 
We develop a novel spatio-temporal infectious disease model called IeRNN, which is a hybrid model consisting of I-equation from SEIR driven by spatial features. With such features and RNN dynamics as external input to the Iequation, the robustness to parameter initialization in model training is greatly improved. In 1-day and 7-day ahead prediction, our model outperforms standard temporal models. In future work, the social control mechanisms (Albi, Pareschi, and Zanella 2020; Morris et al. 2020 ) could be considered to strengthen the I-equation, as well as traffic data to expand inflow effect beyond geographic neighbors. 