Forecasting the COVID-19 transmission in Italy based on the minimum spanning tree of dynamic region network

Theoretical background

The spread of infectious diseases in a region is described as the dynamic evolution of a nonlinear system, while the outbreak of COVID-19 is regarded as a qualitative state transition of a dynamic system (Yang et al., 2020). From the perspective of dynamic modelling, the stage before the COVID-19 outbreak is regarded as a pre-outbreak stage, immediately after which the system undergoes a critical transition. Then the dynamical process of an epidemic system can be roughly modelled as three stages similar to the dynamics of disease progression (Chen et al., 2012): (i) the normal stage, which is a stable stage with high resilience; (ii) the pre-outbreak stage, which is dynamically unstable. At this stage, the epidemic is still controllable through appropriate measures; and (iii) the outbreak stage, which is an uncontrolled stage with high elastic dynamic features. As shown in Fig. 1B, when the system transits from the normal stage to the pre-outbreak stage, the region network changes significantly and the indicator of our model rises sharply. Different from the traditional detection of the outbreak stage, our model could determine the pre-outbreak stage which generally has no obvious abnormalities but with high potential of state transition into a severe and irreversible stage.

Data collection

Data about the COVID-19 pandemic can be obtained from the GitHub repository managed by Johns Hopkins University for modeling (Dong, Du & Gardner, 2020), which contains publicly available data from multiple sources. The open dataset of Italy we use is available on the website COVID-19 in Italy at: https://github.com/pcm-dpc/COVID-19. This database has been created and managed by the Italian Civil Protection Department, which is updated and integrated daily (Italian Civil Protection D et al., 2020).

Algorithm applied in Italy

The MST-DNM model is illustrated in Fig. 2 and the detailed process of our model applied in Italy is described in the following four steps.

i. Modeling and mapping

It is noted that our model is applied to the region network to monitor the COVID-19 transmission and outbreak in Italy. Therefore, it is necessary to construct the regional network based on the Italian regions’ geographic distribution and their adjacent information at first. The adjacent information is shown in Table S1 of Supplementary Information. In the network, each node represents a region or an autonomous province, while each edge represents the adjacent relation between two regions. Then a 21 ∗ 279 data matrix formed by the records of COVID-19 daily confirmed cases is mapped to the network. The region network model of Italy is presented as in Fig. 3.

ii. Weighting and extracting

The region network in the first step can be represented as an undirected graph $G=\left(V,E\right)$, which contains a collection $V={\left\{v_i\right\}}_{i=1}^M$ composed of M vertexes and a collection $E={\left\{e_j\right\}}_{j=1}^N$ composed of N edges in this network. In the Italian region network (Fig. 3), M is 21 and N is 34. It should be noted that in order to make the application of our model easier to understand, Node 15 is not considered in the subsequent calculations. When the model is applied to other regions, the construction of regional network should be reconsidered, instead of simply excluding islands. In addition, as shown in Fig. S1, as the COVID-19 epicenter is located in northern Italy, whether Sardinia is included in the model has little influence on the final warning results.

The number of daily confirmed cases of a region is considered as a sample s to form a set of time series data. Therefore, for each vertex v_i of the region network on the day t, there is a corresponding time series of confirmed cases $S_t^{v_i}=\left\{s_1,s_2,\dots ,s_t\right\}$. In order to assign a weight $W_t^k$ to each edge e^k of the region network on the day t, it is necessary to calculate the correlations between the two vertexes v_i, v_j of the edge e^k as follows:

(1)${\displaystyle W_t^k=\left|\Delta S{D}_t\left(i,j\right)\right|\ast \left|\Delta PC{C}_t\left(i,j\right)\right|,}$

where

(2)${\displaystyle \Delta S{D}_t\left(i,j\right)=\left|\overline{S{D}_t\left(S_t^{v_i}\right)+S{D}_t\left(S_t^{v_j}\right)}-\overline{S{D}_{t-1}\left(S_{t-1}^{v_i}\right)+S{D}_{t-1}\left(S_{t-1}^{v_j}\right)}\right|}$

is the differential standard deviation of the nodes v_i, v_j on day t and t − 1, and $S{D}_t\left(S_t^{v_i}\right)$ and $S{D}_t\left(S_t^{v_j}\right)$ represent the standard deviation of the time series data of the two vertices v_i, v_j. (3)${\displaystyle \Delta PC{C}_t\left(i,j\right)=\left|PC{C}_t\left(S_t^{v_i},S_t^{v_j}\right)\right|-\left|PC{C}_{t-1}\left(S_{t-1}^{v_i},S_{t-1}^{v_j}\right)\right|}$

is the differential Pearson’s correlation coefficient between the two vertices v_i, v_j of the edge e^k, where $PC{C}_t\left(S_t^{v_i},S_t^{v_j}\right)$ and $PC{C}_{t-1}\left(S_{t-1}^{v_i},S_{t-1}^{v_j}\right)$ represent the Pearson’s correlation coefficient between the two vertices v_i, v_j on day t and t − 1 respectively.

Based on the above work, an undirected and edge-weighted network that changes dynamically over time is obtained. The next step is to extract the minimum spanning tree from the dynamic region network at each moment. In detail, when Italy is on the t day of the COVID-19 pandemic, we could extract the minimum spanning tree $MS{T}_t=\left(V,E\right)$ to better describe the evolution of Italian regional network with the change of daily cases. In this work, a classical minimum spanning tree algorithm, the Prim’s algorithm, is applied to the differential weighted network G_t at a specific time t to obtain its minimum spanning tree MST_t.

iii. Calculating early warning indicators I_t

Then, the MST-DNM indicator I_t can be obtained by calculating the weight sum of the minimum spanning tree. $I_t={\sum }_{i=1}^K{W}_t^i$, where $W_t^i$ represents the weight of the edge e_i of the minimum spanning tree MST_t at time t and K is the total number of edges of the minimum spanning tree MST_t. The specific algorithm flow is shown in Algorithm 1.

 
_______________________ 
Algorithm 1 The indicator It using Prim’s algorithm________________________________ 
Require: The nodes of the weighted undirected graph V ; the function w(u,v) 
     which means the weight of the edge (u,v); the function adj(v) which means 
     the nodes adjacent to v. 
Ensure: The sum of weights of the minimum spanning tree in the input graph, 
     It. 
  1:  It ← 0 
  2:  choose an arbitrary node in V to be the root 
  3:  dis(root) ← 0 
  4:  for each node v ∈ (V −{root}) do 
 5:       dis(v) ←∞ 
 6:  end for 
 7:  rest ← V 
 8:  while rest ⁄= Φ do 
 9:       cur ← the node with the minimum dis in rest 
10:       It ← It + dis(cur) 
11:       rest ← rest − cur 
12:       for each node v ∈ adj(cur) do 
13:            dis(v) ← min(dis(v),g(cur,v)) 
14:       end for 
15:  end while 
16:  return It________________________________________________________________________________________    

According to DNB theory, during the critical stage, there are two cases for the minimum spanning tree MST_t at time point t:

• In the MST_t, all of the nodes are DNB members;

• In the MST_t, DNB and non-DNB members both exist.

For the above two cases respectively, the statistical indicator I_t has significant changes as presented in Table 1. Obviously, the MST_t based on the indicator I_t and the edges’ weight W_t has the ability to monitor the dynamical process of COVID-19 spread between regions and issue a warning signal timely.

iv. Identifying early-warning points

In previous studies, machine learning methods, i.e., logistic regression (Yang et al., 2020), have been applied to identify the appearance of critical points based on years of high-dimensional data. However, for COVID-19 which originated at the beginning of 2020, the time series data obtained is of quite small scale, which is difficult for machine learning algorithm to learn the appropriate parameters and features. Therefore, the fold-change threshold, an index of volatility, is used to detect the early-warning signal. Specifically, a 2-fold change threshold is applied to identify the significant changes of the indicator I_t in our study.

The significance of indicator I_t

The transmission of COVID-19 is a complicated dynamic system, which contains many biomedical and social factors. Due to the massive number of influencing factors, it is difficult to describe the transport dynamics in high-dimensional space mathematically. The sharp or qualitative transition of regional network from the normal state to the outbreak state corresponds to the bifurcation point in dynamic system theory (Gilmore, 1993). According to this theory, if the system approaches the critical point, it will eventually be confined to one-dimensional or two-dimensional space, where the dynamic system can be represented in quite simple forms. This is the theoretical basis for developing a general indicator that can describe the dynamics of COVID-19 transmission.

It’s clear from the above statement that the meaning of the variables in Formula (1) is as follows: (i) $\left|\Delta S{D}_t\left(i,j\right)\right|$ can describe the differential fluctuation of cases growth in two adjacent regions compared with the previous time point. (ii) $\left|\Delta PC{C}_t\left(i,j\right)\right|$ can describe the difference of the COVID-19 interaction between two adjacent regions compared with the previous time point. Apparently, attention should be paid to the edge with larger weight. Because it means that the regions associated with this edge not only worsen their own epidemic situation, but also have a great impact on the surrounding regions. Therefore, it’s obvious that the indicator I_t, the weight sum $W_t^k$ of all edges in MST_t, has the ability to observe the change of a group of weighted differential networks.

Early warning of COVID-19 outbreaks in Italy

We collect the historical data of daily cases infected by COVID-19 from February 24, 2020 to November 28, 2020 in Italy. The outbreak points of COVID-19 are defined as the peak of the daily cases.

Provided as in Fig. 4, the early-warning signals are identified through the MST-DNM model for each outbreak of COVID-19. For the first wave of COVID-19 outbreak from mid-February to mid-April, the early-warning signal was issued on March 6, which is about 15 days ahead of the outbreak point. This means that our model has successfully played an early warning role.

On 9 March 2020, the Italian prime minister Mr. G Conte announced the implementation of placing the country into lockdown to restrict the movement of people, thus reducing the possibility of human to human infection (Chintalapudi, Battineni & Amenta, 2020; Remuzzi & Remuzzi, 2020). Since the last week of March, the statistics have become consciously optimistic, and the number of daily cases has been stabilizing. For the second wave of epidemic since early October, later developing into a larger outbreak, the indicator I_t was sensitive and significantly increased about 10 days before the actual number of confirmed cases skyrockets. In addition, the indicator showed a continuous downward trend with wave type after November 5, which means that the number of daily cases in Italy has initially peaked. The successful prediction of each wave of COVID-19 outbreaks in Italy demonstrates the robustness and effectiveness of the MST-DNM model in detecting real-time warning signals for emerging infectious diseases.

The dynamics of COVID-19 transmission in Italy

Dynamic monitoring map

To better illustrate the MST-DNM model’s principle, we introduce the dynamic evolution of the COVID-19 transmission network in Italy. As shown in Fig. 5, the daily number of newly confirmed cases with MinMaxScaler in each region is mapped to each node and the correlation between two vertices of an edge is mapped to the thickness of the edge in the tree network. The specific method of MinMaxScaler is to subtract the minimum value of the feature from the processed value and divide it by the feature range, which is the difference between the original maximum and the original minimum. MinMaxScaler can keep the shape of the original data distribution of each region, and make the coloring of each node unaffected by other regions. It is clear that the edges became thicker before the nodes turned darker on October 20, which indicated our model identified the early-warning signal in the pre-outbreak stage when the actual number of confirmed cases did not increase significantly. After that, i.e., on October 26, the edges continued to become thicker, which meant that the epidemic might continue to worsen in Italy.

Warning function of MST-DNM in Italy

As of November 28, the five regions with the highest cumulative number of confirmed cases in Italy are Lombardy, Piedmont, Campania, Veneto and Emilia-Romagna, which correspond to Nodes 10, 13, 5, 21 and 6 in Fig. 6A, respectively. Among them, Lombardy region is considered as the epicenter of COVID-19 outbreak in Italy (Grasselli, Pesenti & Cecconi, 2020; Tuite et al., 2020).

As shown in Fig. 6A, the dynamic region network is divided into two local networks according to the thickness of the edges, which are centered on Node 10 (Lombardy) and Node 5 (Campania) respectively. It’s obvious that our model has successfully warned two high-risk regions. In fact, a large number of new cases have been confirmed in these regions over the next month as shown in Fig. 6B. In addition, the nodes sandwiched between two local networks, i.e., Node 17 (Tuscany) and Node 19 (Umbria), should also be focused on. It turns out that from October 26 to November 28, the growth rate of new case in Tuscany and Umbria exceeded 200%.

As described in Table 2, the thicker edges in two local networks on October 26, such as $e\left(v_{10},v_{21}\right)$, $e\left(v_6,v_{21}\right)$, $e\left(v_{13},v_{20}\right)$, $e\left(v_5,v_8\right)$ and $e\left(v_{12},v_{14}\right)$, should be focused on. As of November 28, the number of confirmed cases in the regions, like Lombardy, Campania, Lazio and Piedmont, corresponding to the nodes associated with these edges has a high growth rate and a large number of new cases. This also verify the early warning function of MST-DNM, which not only measures its own epidemic situation, but also reflects the regional impact on the surrounding areas.

Table 2:

Description of the early warning function of edge weight.

Local region network #	Edge	Related regions	Newly added cases	Growth rate
1	$e\left(v_{10},v_{21}\right)$	Lombardy	244,726	154.96%
		Veneto	95,506	210.07%
	$e\left(v_6,v_{21}\right)$	Emilia-Romagna	71,379	148.93%
		Veneto	95,506	210.07%
	$e\left(v_{13},v_{20}\right)$	Piedmont	107,150	187.46%
		Valle d’Aosta	3747	140.39%
2	$e\left(v_5,v_8\right)$	Campania	111,077	273.91%
		Lazio	80,262	223.35%
	$e\left(v_{12},v_{14}\right)$	Molise	3193	242.08%
		Puglia	37,341	249.39%

DOI: 10.7717/peerj.11603/table-2

Notes:

“Newly added cases” refers to the cumulative number of confirmed cases in the corresponding region from October 26 to November 28.

“Growth rate” refers to the newly confirmed cases in corresponding regions from October 26 to November 28 divided by the cumulative cases on October 26.

Application of MST-DNM in northern Italy

As shown in Fig. 6B, several areas in northern Italy, such as Lombardy, Veneto, etc., are the most severely affected by COVID-19. Our model can be applied not only to the entire country of Italy, but also to an area. In order to verify the effectiveness of our model, it has also been applied to identify the early-warning signals of COVID-19 outbreaks in northern Italy. The results are presented in Figs. S2 –S3 of Supplementary Information.

Performance comparison

The machine learning algorithms are also used to forecast the COVID-19 epidemic situation (Parhusip, 2020; Singh et al., 2020; Parbat & Chakraborty, 2020). Regarding the identification of early warning signals of COVID-19 outbreaks as a binary classification problem, we compare the performance of our combined model with the support vector machine (SVM). The AUC of MST-DNM is 0.9318, while that of SVM is 0.9076. It’s clear that the performance of a system based on MST-DNM is better than that based on SVM when only the data of daily confirmed cases is provided. In addition, the SVM model issue an early warning signal on September 28, 2020, which is too early to be of practical significance for the second wave outbreak; and our adaptive model can actually issue an early warning signal about 10 days before the actual number of confirmed cases skyrockets. Actually, compared with traditional machine learning algorithms, the MST-DNM model has the following internal strengths. First of all, it is a model-free approach without any training and testing processes. There is no feature selection in MST-DNM strategy, which solely depends on the statistical conditions of our model. Second, it’s noted that there is no limitation of the data sample size for our approach, which means that our model could achieve a good performance with only small sample data. Therefore, it can be applied to describe and monitor the emerging infectious diseases like COVID-19. In addition, our combined model is capable to describe the dynamic process of the spread of COVID-19 through the minimum spanning tree of dynamic region networks.