Real-Time Differential Epidemic Analysis and Prediction for COVID-19 Pandemic

The novel coronavirus (COVID-19) epidemic is generating significant social, economic, and health impacts and has highlighted the importance of real-time analysis and prediction of emerging infectious diseases and health care resource and personnel planning and economical activity guideling.

One of the well-known models that is reasonably predictive for human-to-human transmission is the so-called Susceptible-Infectious-Removed (SIR) model, which was published in its first form around the 1920s [21]. The model later has been extended to consider more complicated situations such as the Susceptible-Exposed-Infected-Removed (SEIR) model [3]. Those models are very successful to describe how the disease dynamics will change over time once more sufficient data are available (from outbreak to finish). Recently those models have been applied to study recent COVID-19 transmissions in different countries with certain successes [16, 4, 17, 25, 15]. For example, it is observed that city-wide lockdown can lower the transmission rate substantially from those models. On the other hand, the data-driven and curve-fitting methods for the prediction of COVID-19 such as Gaussian function based fitting method in IHME projection [10], exponential curving fitting in[24], machine-learning based approaches in [23, 9] can fit the data well. However, those methods generally suffer the lack of physical insights of transmission and will not work well when the data does fit their models well (For instance, two waves of infections). For real time prediction, large projection ranges have to be given, which renders the projection less valuable for medical resources and mitigation policy planning.

Furthermore, the worldwide public health crisis like 1918 pandemic, which killed an estimated 50 million people worldwide, including an estimated 675,000 people in the United States, are difficult to model, even to project [2]. H1N1 flu lasted from April 1917 to April 1919 for two years with three major phases across different parts of world. Cities like San Francisco even experienced strong second waves of death increase when the social distance was relaxed too early. Such long-term multi-phase and multi-year transmission dynamic is very difficult to be captured by existing pandemic models. As a result, real-time short-term prediction and near-tern trend projection can give each city and county policy maker and resource planners extremely valuable information to guide the disease mitigation and business open/close decision in a real time.

On the other hand, traditional physics-based SIR/SEIR epidemic modeling and its variants suffer several drawbacks especially for real-time disease transmission predictions as it is static model in which many key parameters such as transmission rate and recovery rates are typically fixed values from fitting. However, for countries like US, the local government interventions and mitigation policies such as active surveillance, contact tracing, quarantine, massive testing, school and business closure, shelter in place, social distancing etc. keep changing to respond the local epidemic situations. Also for each state in US has different time lines for implementing different prevention and mitigation measures, which make the static based prediction even more difficult to predict. Further the social behaviors of population, such as bearing face coverage masks are not consistent through different regions and cities as the pandemic progresses, which will affect the transmission rate as well. The availability of healthcare resources of each city and county, which are also changing factors, which may also affect recovery rate. As a result, existing static SEIR models do not fit well for real-time disease prediction as the key parameters such as transmission rate, recovery rate, which is also closely related to the effective reproduction rate $R0$, are time-varying parameters. On the other hand, one important observation is that contagious disease transmission prediction is similar to weather prediction: only short term prediction is accurate due to many time-varying interplaying factors as well as social and behavior uncertainties involved. The widely watched IHME (Institute for Health Metrics and Evaluation) model for real-time COVID-19 epidemic prediction [10] is purely based on mathematic curve fitting techniques, which lacks the theoretical foundation of epidemic transmission found in the SIR models and thus its prediction accuracy is highly debatable.

To mitigate this problem, time-varying SEIR models have been proposed in the past. Dureau et al tried to model the time-varying affects of the SEIR models by considering partial and noisy data. They introduced stochastic processes into the SEIR models and solved the resulting stochastic SEIR partial differential equations using Markov Chain Monte Carlo methods in which the transmission is modeled as random walks, which are very expensive to compute [8]. This approach was applied to study the transmission within and outside Wuhan for January to February 2020 [11]. Recently Chen et al proposed time-dependent SIR model to model the COVID-19 outbreaks in China [6], In this mode, model parameters are computed in a daily basis. As a result, it lacks the good predictability as all the key parameters have been predicted a prior first, and some ad-hoc methods was introduced to predict the parameters in a near future.

In this work, we propose a new real-time differential virus transmission modeling method, which can give more accurate and robust short-term predictions of COVID-19 transmitted infectious disease while still maintain the near-term projection benefit. The new model is based on enhanced Susceptible-Exposed-Infected-Removed (SEIR) virus transmission model. But it tries to obtain the differential view of pandemic dynamics in a short history window to analyze the short trend of the transmission so that it can be more accurate over different periods of time to accommodate the adaptive changes of disease mitigation, business activity and social behavior of populations. As the parameters of the improved SEIR models are trained by short history window data for accurate trend prediction, our differential epidemic model, essentially are window-based time-varying SEIR model. Since SEIR model still is a physics-based disease transmission model, its near-term (like one month) projection can still be very instrumental for policy makers to guide their decision for disease mitigation and business activity policy changes in a real-time. Numerical results on the recent COVID-19 data from China, Italy and US, California and New York states are analyzed. A dedicated website has been built to show the projections based on the latest data [1].

Fig. 1 illustrates the modeling and prediction results from the proposed differential SEIR model for the COVID-19 disease for California from early March to later July 2020. As we can see, at the beginning, the infected grows very fast, the projected growth also reflects such fast changing rates. But as social distancing and stay at home policies were introduced across cities and counties in California in the later March, the growth rate went down, the projected growth at different time point also reflects such trends. Based on the projection of our model at May 1, 2020, the infected case will reach to the peak about 72.37K around Jun 14 in California, which indicates the peak medical resources needed. In contrast, the well-watched IHME’s prediction [10] predicts that the peak medical resource needed is around April 17. We will show more significant differences between our models and IHME’s prediction  [10] for New York state later.

In this section, we first present the extended SIR model, which is a extension of classical SEIR model [20, 19, 12, 18, 7, 22, 13] as the base model for the proposed differential modeling shown later. The proposed base model, called SEIRDP model, is similar to the recently proposed generalized SEIR model for studying the COVID-19 disease in Wuhan and China [17]. We removed the quarantine compartment in the proposed SEIRDP model as there is no quarantine data for most of countries outside China. The resulting SEIRDP (Susceptible, Exposed, Infected, Recovered, Death, Insusceptible (P)) model is shown in Fig. 2. In this base model, we have six states, i.e. $\{S(t),P(t),E(t),I(t),R(t),D(t)\}$, which indicate the number of the susceptible cases, insusceptible cases, exposed cases (infected but not yet be infectious, in a latent period), infectious cases (confirmed with infectious capacity), recovered and immune cases and closed cases (or death). Then the total number of population in a certain region or county is $N=S+P+E+I+R+D$. The coefficients $\{\alpha,\beta,\gamma^{-1},\lambda^{-1},\kappa\}$ represent the protection rate, infection rate, average latent time, cure rate, and mortality rate, respectively. The introduction of insusceptible compartment represents gradual changing (or growing) population, which will not be infected due to some strong disease mitigation measures such as enforced shelter in place, strict city lockdown in China and European countries. The basic reproduction number, $R0$, represents the number of secondary infections from a primary infected individual in a fully susceptible population, which can be computed by [17]

$$R0=\beta\lambda^{-1}(1-\alpha)^{n}$$

(1)

where $n$ is the number of days. We remark that If we force $\alpha=0$, the proposed SEIRDP model essentially becomes the classical SEIR model (if we put recovered and death cases into one recovery compartment).

Once we have the base SEIR model, then we can present our differential SEIR model. In this model, basically the key parameters will become time dependent.

In our work, we introduce a time window concept, which represents a short period of time in a few days, as shown in Fig. 3. In this figure, the red box present a history length of the data ($h$) we use for training the model. The red box will move one window size ($w$) as one step forward. For instance five days can be a good time window. In each time window, all the key parameters in the base SEIR models are constant. This reflects the fact that the transmission and recovery conditions for population does not change in a short window time so that we can use a traditional static SEIR model based on the history of data around this window time and we can also predict a short future with sufficient accurate assuming such disease dynamic does not change dramatically, which is also reasonable. Based on this observation, we can present the resulting ordinary differential equations for the proposed differential SEIR model as follows:

$$\begin{split}&\frac{dS(t)}{dt}=-\beta(t_{i})\frac{S(t)I(t)}{N}-\alpha(t_{i})S(t)\\ &\frac{dE(t)}{dt}=\beta(t_{i})\frac{S(t)I(t)}{N}-\gamma(t_{i})E(t)\\ &\frac{dI(t)}{dt}=\gamma(t_{i})E(t)-\lambda(t_{i},t)I(t)-\kappa(t_{i},t)I(t)\\ &\frac{dR(t)}{dt}=\lambda(t_{i},t)I(t)\\ &\frac{dD(t)}{dt}=\kappa(t_{i},t)I(t)\\ &\frac{dP(t)}{dt}=\alpha(t_{i})S(t)\end{split}$$

(2)

where time $t$ belongs to the $i$th time window. If we use $t_{i}$ to indicate the last day in the window $i$, then we have $t\in[t_{i}-h+1,t_{i}]$, where $h$ is the window size. The window can be moved forward with stride $w$ days. Therefore, we have the relationship $t_{i+1}=t_{i}+w$. The cure rate $\lambda(t_{i},t)$ and mortality rate $\kappa(t_{i},t)$ are time-dependent, which again are dependent on two parameters. They again are time-window dependent (explained below):

$$\begin{split}&\lambda(t_{i},t)=\lambda_{0}(t_{i})[1-\exp(-\lambda_{1}(t_{i})t)]\\ &\kappa(t_{t},t))=\kappa_{0}(t_{i})\exp(-\kappa_{1}(t_{i})t)\\ \end{split}$$

Those two equations basically say that cute rate or recovery rate $\lambda(t)$ will goes to constant value $\lambda_{0}$ exponentially over a time and the mortality rate $\kappa(t)$ goes down exponentially over time [5].

We notice that all the key parameters for the enhanced SEIR model are time-varying instead of fixed values. But they are so-called time window dependent as they do not change every day compared to existing time-varying based SEIR models. Specifically, we select $w=5$ days as a window (the window size is a hyper parameter and can be optimized for different regions and countries). We perform each prediction over each non-overlapping window and then slide one window forward for next prediction and so on. As a result, within the $i$th window, the seven parameters, $\{\alpha(t_{i})$, $\beta(t_{i})$, $\gamma^{-1}(t_{i})$, $\lambda_{0}^{-1}(t_{i})$, $\lambda_{1}^{-1}(t_{i})$, $\kappa_{0}(t_{i})$,$\kappa_{1}(t_{i})\}$ are kept constant. The parameters will be found by a regression process shown next section. The stride can be $h=10$ days or $h=15$ days, which is also another hyper parameter for our model. Since all the parameters are time window dependent, the reproduction number becomes time dependent, which is also called effective reproduction number, $R_{t}$ at time $t$ in $i$th time window. $R_{t}$ can be computed as follows:

$$R_{t}(t_{i})=\beta(t_{i})\lambda^{-1}(t_{i})\prod_{j=1}^{n}(1-\alpha_{j})$$

(3)

We note that such window-based SEIR model obtained at $i$th window also depends on the historical data from the previous windows. The reason is that our differential model can be viewed as performing the differential operations on a dynamic systems at a specific time frame or window, not on a static function. As a result, the impacts from the historical data of the previous windows will be represented as the initial conditions for solving the PDE of (2) and obtaining the resulting parameters by fitting the resulting discretized PDE with the data in this window.

For the proposed differential SEIR model, for each differential window, we will calculate the seven $\{\alpha(t_{i})$, $\beta(t_{i})$, $\gamma^{-1}(t_{i})$, $\lambda_{0}^{-1}(t_{i})$, $\lambda_{1}^{-1}(t_{i})$, $\kappa_{0}(t_{i})$,$\kappa_{1}(t_{i})\}$ in the $i$th time window indicated by $t_{i}$, which indicates the last day in the window $i$. As a result, we can rewrite the partial differential equation (2) into the following initial value PDE in matrix form with the initial condition $\mathbf{Y}_{0}$, which indicate the the impacts from the historical data before $i$th window:

$$\begin{split}&\mbox{PDE}:\frac{d\mathbf{Y}}{dt}=\mathbf{A}\cdot\mathbf{Y}+\mathbf{F}\\ &\mbox{IC}:\mathbf{Y}_{0}\end{split}$$

(4)

where IC is initial condition.

$$\begin{split}&\mathbf{Y}=\begin{bmatrix}S(t)&E(t)&I(t)&R(t)&D(t)&P(t)\\ \end{bmatrix}^{T}\\ &\mathbf{A}=\begin{bmatrix}-\alpha(t_{i})&0&0&0&0&0\\ 0&-\gamma(t_{i})&0&0&0&0\\ 0&\gamma(t_{i})&-\lambda(t_{i},t)-\kappa(t_{i},t)&0&0&0\\ 0&0&\lambda(t_{i},t)&0&0&0\\ 0&0&\kappa(t_{i},t)&0&0&0\\ \alpha(t_{i})&0&0&0&0&0\\ \end{bmatrix}\\ &\mathbf{F}=S(t)\cdot I(t)\cdot\begin{bmatrix}-\frac{\beta(t_{i})}{N}&\frac{\beta(t_{i})}{N}&0&0&0&0\\ \end{bmatrix}^{T}\\ &\mathbf{Y}_{0}=\mathbf{Y}(t=t_{i}-h+1)=\mathbf{Y}(t=t_{i-1}-h+w+1)\end{split}$$

Then we perform the time domain discretization using simple Forward Euler or higher order explicit Runge-Kutta method and we end up of number of algebraic equations. By solving this resulting equation in the time domain, we can obtain the seven parameters over the given data of $\{I(t),R(t),D(t)\}$ in the past 10 or 15 days with the initial state conditions computed from the previous time window. In this paper, we use a nonlinear least-squares solver to estimate the parameters by the expression [5]

$$\hat{\theta}=\arg\min_{\theta}:=||\text{SEIRDP}(\theta,time)-\begin{bmatrix}I\\ R\\ D\\ \end{bmatrix}||^{2}$$

(5)

where $\theta=\{\alpha(t_{i})$, $\beta(t_{i})$, $\gamma^{-1}(t_{i})$, $\lambda_{0}^{-1}(t_{i})$, $\lambda_{1}^{-1}(t_{i})$, $\kappa_{0}(t_{i})$,$\kappa_{1}(t_{i})\}$, $\hat{\theta}$ is the parameters of estimated model, $time$ is the time from the real data, $\{I,R,D\}$ are the infected, recovered and death cases from the real data, and SEIRDP$(\cdot)$ represents the SEIRDP model.

In this paper, we have proposed a new real-time differential virus transmission model, which can give more accurate and robust short-term predictions of COVID-19 transmitted infectious disease with benefits for near-term trend projection. The new model is based on enhanced Susceptible-Exposed-Infected-Removed (SEIR) virus transmission model. As the parameters of the improved SEIR models are trained by short history window data for accurate trend prediction, our differential epidemic model, essentially are window-based time-varying SEIR model. Numerical results on the recent COVID-19 data from China, Italy and US, California and New York states have been analyzed.