key: cord-0465155-5uvv46mi
authors: Prasse, Bastian; Achterberg, Massimo A.; Ma, Long; Mieghem, Piet Van
title: Network-Based Prediction of the 2019-nCoV Epidemic Outbreak in the Chinese Province Hubei
date: 2020-02-11
journal: nan
DOI: nan
sha: 70565fd845b58e8d897aed13c8d855c82bd21e5f
doc_id: 465155
cord_uid: 5uvv46mi

At the moment of writing (12 February, 2020), the future evolution of the 2019-nCoV virus is unclear. Predictions of the further course of the epidemic are decisive to deploy targeted disease control measures. We consider a network-based model to describe the 2019-nCoV epidemic in the Hubei province. The network is composed of the cities in Hubei and their interactions (e.g., traffic flow). However, the precise interactions between cities is unknown and must be inferred from observing the epidemic. We propose a network-based method to predict the future prevalence of the 2019-nCoV virus in every city. Our results indicate that network-based modelling is beneficial for an accurate forecast of the epidemic outbreak.

In December 2019, the novel coronavirus 2019-nCoV emerged in the Chinese city Wuhan [1] . Individuals that are infected by the 2019-nCoV virus suffer from the Novel Coronavirus Pneumonia (NCP). Contrary to initial observations [2] , the 2019-nCoV virus does spread from person to person as confirmed in [3] . On February 12, 2020, there were more than 45,000 confirmed infections, and more than 1000 people died [4, 5, 6] . Assessing the further spread of the 2019-nCoV epidemic poses a major public health concern.

Many studies aim to estimate the basic reproduction number R 0 of the 2019-nCoV epidemic [7, 8, 9, 10, 11, 12, 13, 14, 15] . The basic reproduction number R 0 is a crucial quantity to evaluate the hostility of a virus [16, 17] . The basic reproduction number R 0 is defined [18] as "The expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual during its entire period of infectiousness". The greater the basic reproduction R 0 , the more individuals are infected in the long-term endemic state of the virus. If R 0 < 1, then the virus dies out. The estimates for the basic reproduction number R 0 of the 2019-nCoV epidemic range from R 0 = 2.0 to R 0 = 3.77.

The basic reproduction number R 0 only coarsely assesses the quantitative behaviour of the epidemic. To obtain a more detailed picture of the epidemic, the development of epidemic outbreak prediction methods is focal. A diverse body of research considers the prediction of general epidemics. For instance, prediction methods are based on Kalman filtering [19] , Bayesian model averaging [20] , basic regression [21] and kernel density estimation [22] . Recent work focussed on the dependency of population flow and the viral spread [23, 24, 25, 26] . As shown by Pei et al. [27] , the spread of influenza can be more accurately predicted by taking the population flow between cities into account. Read et al. [14] predicted the 2019-nCoV epidemic by using the Official Aviation Guide (OAG) Traffic Analyser dataset. Additionally to the OAG dataset, Wu et al. [15] used the Tencent database to predict the 2019-nCoV viral spread.

The population flow clearly has an impact on the evolution of an epidemic. However, the exact population flow is unknown, and epidemic prediction methods must account for inaccuracies of population flow data. In this work, we consider the most extreme case by assuming no prior knowledge of the population flow. To forecast the 2019-nCoV epidemic, we design a network-based prediction method that estimates the interactions between cities as an intermediate step. On February 9th, 2020, approximately 70% of the global 2019-nCoV infections are located in the Chinese province Hubei. Thus, we focus on the 2019-nCoV epidemic in Hubei. Our goal is to predict the 2019-nCoV outbreak for every city in Hubei. Section 2 introduces the available data on the 2019-nCoV virus in Hubei. The epidemic model is proposed in Section 3, and the prediction method is outlined in Section 4. The prediction accuracy is evaluated on past data in Section 5.

The time series of reported infections in Hubei forms the basis for the epidemic outbreak prediction. Hubei is divided into 17 cities (more precisely, prefecture-level divisions) and contains the city Wuhan, as illustrated by Figure 1 . We do not consider the city Shennongjia, since the number of infections in Shennongjia is small. We denote the number of considered cities by N = 16. The number of newly reported infections for each city in Hubei is openly accessible via the website of the Hubei Province Health Committee [28] . The data is updated daily and follows the standard time offset of UTC+08:00. Except for Wuhan, the total number of reported infections is small before January 21, 2020. Hence, we consider the 2019-nCoV epidemic outbreak starting from January 21. We denote the discrete time by k ∈ N. The difference of time k to k + 1 equals one day, and the initial time k = 1 corresponds to January 21, 2020. The website [28] states the number of reported infections N rep,i [k] at every time k in every city i = 1, ..., N . We obtain the population size p i of each city i from the Hubei Statistical Yearbook [29] . The reported fraction of infected individuals in city i at time k follows as 

We model the spread of the 2019-nCoV virus by the SIR-model: At any discrete time k, every individual is in either one of the compartments susceptible (healthy), infectious or removed. Susceptible individuals can get infectious due to contact with infectious individuals. Due to hospitalisation, quarantine measures or death, infectious individuals become removed individuals, which cannot infect 

(1)

and the fraction of susceptible individuals follows as

Here, β ij denotes the infection probability from city j to city i, and δ i denotes the curing probability of city i.

The SIR model (1) assumes that the spreading parameters δ i , β ij do not change over time k. The curing probability δ i quantifies the capacity of individuals in city i to cure from the virus. The infection probability β ij specifies the number of contacts of individuals in city j with individuals in city i. We emphasise that β ii = 0 since individuals within one city i do interact with each other. The contact network between cities in Hubei is given by the N × N matrix

whose elements are probabilities 0 ≤ β ij ≤ 1. Neither the curing probabilities δ i nor the infection probabilities β ij are known for the 2019-nCoV epidemic. Potentially, it is possible to state bounds or estimates for the spreading parameters δ i and β ij by making use of the people flow or geographical distances between the respective cities. Nevertheless, there would remain an uncertainty regarding the precise value of the spreading parameters δ i and β ij . In this work, we consider the most extreme case: there is no a priori knowledge on the curing probabilities δ i nor the infection probabilities β ij .

In Section 4, we develop an inference method to estimate the spreading parameters δ i and β ij from observing the epidemic.

We propose a network-based method to predict the outbreak of 2019-nCoV virus, which consists of three steps. First, we preprocess the raw data of the confirmed number of infected individuals in Subsection 4.1 to obtain an SIR time series v i [1] , v i [2] , ... of the viral state for every city i. Second, based on the time series v i [1] , v i [2] , ..., we obtain estimatesδ i andβ ij of the unknown spreading spreading parameters δ i and β ij in Subsection 4.2. Third, the estimatesδ i andβ ij result in an SIR model (1), which we iterate for future times k to predict the evolution of the 2019-Cov virus. Subsection 4.1 and Subsection 4.2 give an outline of the first two steps of the prediction method. We refer the reader to Appendix B for a detailed description of the prediction method.

We denote the number of observations by n, which equals the number of days since January 21, 2020. Our goal is to obtain an SIR viral state vector

) T for every city i at any time k = 1, ..., n based on the data described in Section 2. The fraction of susceptible individuals follows as (1), if the curing probability δ i were known. However, we do not know the curing probability δ i . Hence, we consider 50 equidistant candidate values for the curing probability δ i , ranging from δ min = 0.01 to δ max = 1. We define the set of candidate values as Ω = {δ min , ..., δ max }. For every candidate value δ i ∈ Ω, the fraction of removed individuals R i [k] follows from (1) at all times k ≥ 2. Thus, we obtain 50 potential sequences R i [1] , ..., R i [n], each of which corresponding to one candidate value δ i ∈ Ω. We estimate the curing probability δ i , and hence implicitly the sequence

, as the element in Ω that resulted in the best fit of the SIR model (1) 

For every city i, the curing probability δ i is estimated as one of the candidate values in Ω, as outlined in Subsection 4.1. The remaining task is to estimate the infection probabilities β ij . The goal of network inference [33, 34, 35, 36] is to estimate the matrix of infection probabilities B from the SIR viral state

. The matrix B can be interpreted as a weighted adjacency matrix. We adapt a network inference approach 4 [38, 31] , which is based on formulating a set of linear equations and the least absolute shrinkage and selection operator (LASSO) [39, 40] . The crucial observation from the SIR governing equations (1) is that β ij appears linearly, whereas the state variables S i , I i and R i do not. From (1), the infection probabilities β ij satisfy

for all cities i = 1, ..., N . Here, the (n − 1) × 1 vector V i and the (n − 1) × N matrix F i are given by

. . .

3 Potentially, the outlier is due to the increase in the maximum number of individuals that can be diagnosed in Wuhan, from 200 to 2000 individuals per day as of January 27th [32]. 4 The network inference approach [31] is also applicable to general compartmental epidemic models [37] , such as the Susceptible-Exposed-Infected-Removed (SEIR) epidemic model. and

If the SIR model (1) were an exact description of the evolution of the coronavirus, then the linear system (2) (2) only holds approximately. Thus, we resort to estimating the infection probabilities β ij by minimising the deviation of the left side and the right side of (2). We reconstruct the network by the LASSO [39, 40] as follows:

The first term in the objective function of (5) measures the deviation of the left side and the right side of (2). The sum in the objective of (5) is an 1 -norm regularisation term which avoids overfitting. We choose to not penalise the self-infection probability β ii , since we expect the infections among individuals within the same city i to be dominant. The regularisation parameter ρ i > 0 is set by cross-validation. The LASSO network inference (5) allows for the incorporation of a priori knowledge of the contact network B by adding further constraints to the infection probabilities β ij . We emphasise that an accurate prediction of an SIR epidemic outbreak does not require an accurate network inference [31] .

The accuracy of the network-based prediction method in Section 4 is evaluated by comparison to a simple prediction method. Qualitatively, the virus spread in many epidemiological model follows a sigmoid function, see also [42] . A particular sigmoid function is obtained by logistic regression. As a comparison to the method in Section 4, we apply logistic regression on the reported fractions I rep,i [1] , ..., I rep,i [n] of infection individuals, independently for each city i in Hubei. Logistic regression is advantageous because a logistic function is a closed-form expression, and its parameters can be determined by non-linear regression. Moreover, the logistic function is an approximation to the exact solution of some epidemiological models and population growth models [41, 42, 43] . For further details regarding logistic regression, we refer the reader to Appendix C. We denote the cumulative fraction of infections at time k by

At the time of writing, the data is available from January 21 until February 11, 2020. To evaluate the prediction accuracy, we remove the data for a fixed number of days, say m, prior to February 11. The prediction model is determined upon the data from 21 January up to 11 − m February, 2020. Then, we predict the course of the disease up to February 11, and the number of omitted days m is equal to the number of prediction days. The course of the disease is shown in Figure 2 For most predictions shown in Figure 2 , the logistic curve appears to underestimate the true fraction of infected individuals, whereas the network-based method seems to overestimate the true value. The logistic curve is therefore a lower bound prediction for the real fraction of infected individuals.

The prediction accuracy decreases if the prediction time is increased, which we quantify by the Mean Absolute Percentage Error (MAPE)

at any time k. Here,Î cs,i [k] denotes the predicted cumulative fraction of individuals of city i at time k. Figure 3 depicts the MAPE prediction error for the data shown in Figure 2 . Two observations are worth mentioning. First, as expected, the prediction error increases when predicting more days ahead. Second, the network-based method always provides more accurate predictions than the logistic regression. Figure 4 illustrates the prediction accuracy versus the time that the epidemic outbreak has been observed. As the epidemic evolves, the prediction accuracy increases.

Finally, we consider the prediction of the fraction of infected individuals for the next five days. We stress that, as shown in Figure 2 and Figure 4 , the prediction might be inaccurate for more than four days ahead. The predicted number of infected individuals for each city i is shown in Table 1 .

We applied a network-based SIR epidemic model to predict the outbreak of the 2019-nCoV virus for each city in the Chinese province Hubei. The epidemic model allows to explicitly specify the interactions of individuals of different cities, for instance by using traffic patterns between cities. However, the precise interactions between cities is unknown and must be inferred from observing the evolution of the epidemic.

We proposed a network-based prediction method, which estimates the interactions between cities as an intermediate step. We did not assume any prior knowledge on the interactions between cities. The prediction method is evaluated on past data of the 2019-nCoV outbreak in Hubei. Our results indicate that a network-based modelling approach may yield more accurate predictions than modelling the epidemic for each city independently. We believe that the prediction method can be further improved, e.g., by using traffic flow patterns as prior knowledge. True Data Predicted (network-based) Predicted (logistic regression) Figure 5 : The predicted number of infected individuals for 4 cities in Hubei by both methods for the coming 5 days. The first prediction data point (February 11, 2020) coincides with the last day that has been observed. We selected 4 of the N = 16 cities to avoid overlapping curves. The predicted number of infected individuals for the network inference method is presented in Table 1 . A Details on the Data of the 2019-nCoV Epidemic Outbreak Table 2 shows the cities of the province Hubei and the respective population size p i for every city i.

The time series of the reported number of infections N rep,i [k] is stated in Table 3 . for the infection probabilities β i1 (δ i ), ..., β iN (δ i ). Furthermore, the network inference returns the mean squared error MSE (δ i ), which corresponds to the first term in the objective of (5). The smaller the mean squared error MSE (δ i ), the better the fit of the SIR model (1) To determine the regularisation parameter ρ i in the LASSO (5), we consider 100 candidate values specified by the set Θ i = {ρ min,i , ..., ρ max,i }. In line 4 of Algorithm 2, the maximum value is set to [44] the solution to the LASSO (5) is β ij = 0 for all cities j. For every value of the regularisation parameter ρ i ∈ Θ i , we compute the mean squared error MSE (δ i , ρ i ) by 3-fold cross-validation [40] . For every fold, the rows of the matrix F i and the vector V i are divided into a training set F i,train , V i,train and a validation set F i,val , V i,val . We compute the solution β i1 , ..., β iN to the LASSO (5) on the training set of every fold F i,train , V i,train . The mean squared error MSE (δ i , ρ i ) then equals

averaged over all folds. Finally, we set the regularisation parameter ρ i to the minimiser of MSE (δ i , ρ i ). 

A novel coronavirus emerging in China-key questions for impact assessment

A familial cluster of pneumonia associated with the 2019 novel coronavirus indicating person-to-person transmission: a study of a family cluster

Pattern of early human-to-human transmission of Wuhan 2019-nCoV

Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

The mathematics of infectious diseases

A brief history of R 0 and a recipe for its calculation

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Comparison of filtering methods for the modeling and retrospective forecasting of influenza epidemics

Individual versus superensemble forecasts of seasonal influenza outbreaks in the United States

epiforecast: Tools for forecasting semi-regular seasonal epidemic curves and similar time series

Prediction of infectious disease epidemics via weighted density ensembles

The role of the airline transportation network in the prediction and predictability of global epidemics

Multiscale mobility networks and the spatial spreading of infectious diseases

Natural human mobility patterns and spatial spread of infectious diseases

The hidden geometry of complex, network-driven contagion phenomena

Forecasting the spatial transmission of influenza in the United States

Hubei Statistical Yearbook

An individual-based approach to SIR epidemics in contact networks

Network reconstruction and prediction of epidemic outbreaks for general group-based compartmental epidemic models

Inferring network properties based on the epidemic prevalence

Network inference from population-level observation of epidemics

Revealing networks from dynamics: an introduction

Data based identification and prediction of nonlinear and complex dynamical systems

Generalized epidemic mean-field model for spreading processes over multilayer complex networks

Network reconstruction and prediction of epidemic outbreaks for NIMFA processes

Regression shrinkage and selection via the lasso

Statistical learning with sparsity: the lasso and generalizations

Notice sur la loi que la population suit dans son accroissement

Universality of the SIS prevalence in networks

Time-dependent solution of the NIMFA equations around the epidemic threshold

An interior-point method for largescale l 1 -regularized least squares

Algorithm 1 Epidemic outbreak prediction 1: Input: reported fraction of infections I rep,i [1], ..., I rep,i [n] for all cities i; prediction time n pred 2: Output: predicted fraction of infectionsÎ i

← smoothdata(I rep,i [1], ..., I rep,i [n]) for all i = 1

δ i ) , ..., β iN (δ i ) , MSE (δ i )) ← Network inference(δ i

,β iN ) ← β i1 (δ i ), ..., β iN (δ i ) 16: end for Step 3 -Iterating SIR model 17: for i = 1

β iN (δ i ) for the infection probabilities is obtained by solving the LASSO (5) on the whole matrix F i and vector V i . To solve the LASSO (5) numerically

We are grateful to Fenghua Wang for helping with collecting the data. Long Ma is grateful for the support from the China Scholarship Council.

estimate MSE(δ i , ρ i ) by 3-fold cross-validation on F i , V i and solving (5) on the respective training set 9: end for 10: ρ opt,i ← argmin

A logistic curve is given by the following equationIn our formulation, y(t) is the time-dependent fraction of infectious individuals, t is the time in days, where January 21 serves as initial condition (t = 0), y ∞ is the fraction of infected individuals when time approaches infinity, K is the logistic growth rate and t 0 indicates the inflection point of the logistic equation. For each city in Hubei, we have applied the Matlab command lsqcurvefit to fit the reported cumulative fractionof infected individual to equation (6).