Spatio-temporal distribution pattern of COVID-19 in the Northern Italy during the first-wave scenario: The role of the highway network.

Background: The rapid outbreak of Coronavirus disease 2019 (COVID-19) has posed several challenges to the scientific community. The goal of this paper is to investigate the spread of COVID-19 in Northern Italy during the so-called first wave scenario and to provide a qualitative comparison with the local highway net.
Methods: Fixed a grid of days from February 27, 2020, the cumulative numbers of infections in each considered province have been compared to sequences of thresholds. As a consequence, a time-evolving classification of the state of danger in terms of Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infections, in view of the smallest threshold overtaken by this comparison, has been obtained for each considered province. The provinces with a significant amount of cases have then been collected into matrices containing only the ones featuring a significant amount of cases.
Results: The time evolution of the classification has then been qualitatively compared to the highway network, to identify similarities and thus linking the rapid spreading of COVID-19 and the highway connections. Conclusions: The obtained results demonstrate how the proposed model properly fits with the spread of COVID-19 along with the Italian highway transport network and could be implemented to analyze qualitatively other disease transmissions in different contexts and time periods.

Introduction

Background

The aim of this paper is to study the spread of Coronavirus Disease-19 (COVID-19) in the Northern Italy during the so-called first wave scenario. For the sake of the clarity, it is important to recall that Italy is divided by the administrative point of view into 20 regions, each of them considered as an autonomous entity with defined powers, including the financial control of the healthcare system (Torri et al., 2020). They can be grouped following their geographical distribution, into regions of North, Center and South of Italy. With the only exception of Val d’Aosta, each region is divided into a various number of second-level administrative entities, named provinces, which can be viewed as an intermediated level between regions and municipalities and that, in general, are named after their biggest or most representative cities, the so-called province capital cities.

In late February 2020, a sudden surge in Coronavirus Disease-19 (COVID-19) occurred in the Northern Italy. When the first Italian cases have been detected on February 21, 2020, in the municipality of Codogno, belonging to the Lodi province in the Lombardia region, and in Vò Euganeo, Padova province in the Veneto region (Fig. 1 ) (Gianquintieri et al., 2020). The total number of people confirmed to be infected by the coronavirus SARS-CoV-2 was 16 (Indolfi and Spaccarotella, 2020).

Fig. 1

First outbreacks in Italy (Lab24, Il Sole 24 ore).

On February 23, the number of infections leapt to 152, a jump of around 100 within a day (Italian Civil Protection, 2020. In the following days, a massive increase in the number of new infections has been observed in the northern Italy. From February 24 on, the surfaced outbreaks lead to lockdown measures for 10 municipalities in the province of Lodi and one in the province of Padova (Bertuzzo et al., 2020). On March 8, other containment measures were announced for the whole country to counterbalance the Italy’s severe COVID-19 outbreak with the purpose to obtain an effective reduction in transmission after a time lag of 14–18 days (Vinceti et al., 2020). More in detail, strict restrictions regarding travels and individual mobility of citizens have been adopted (Chinazzi et al., 2019, Sebastiani et al., 2020), since analogous measures were already applied successfully in China (Kraemer et al., 2020) since January 2020.

Having regard to the Italian scenario, several attempts have been made to study the connection between the spread of the SARS-CoV-2 (from now on, coronavirus or virus) and various transport systems (Cartenì et al., 2021, Coelho et al., 2020, Gondauri and Batiashvili, 2020).

The hypothesis that the diffusion of the coronavirus during its early phase in Lombardia has followed the Italian highway network has been formulated by Sebastiani (2020) (see also Usuelli, 2020, Micheli et al., 2020). Similarly, Maietti et al. (2020) analyzed the COVID-19 infection rate per 1000 inhabitants for the nine provinces of the Emilia-Romagna region on May 25, suggesting that the ancient Roman road Via Emilia (also named as SS9) played a pivotal role in the spread of the virus throughout Emilia Romagna. Giuliani et al. (2020) proposed an endemic-epidemic time-series mixed-effects generalized linear model for areal disease counts: they found that containment policies limited the spread to nearby areas. D’Angelo et al. (2021) investigated the spatio-temporal spread pattern of COVID-19 in Italy from February to October 2020: they confirmed the conclusions of Giuliani et al. (2020), and suggested that suggests the temporal evolution of cases in each province does not depend on the temporal evolution of the other provinces.

Moreover, Gatto et al., 2020, Amer and Bergquist, 2021 proposed a spatially explicit epidemiological model, mentioning the link between the radiation of the epidemic along transport network. Paez et al. (2021) carried out a spatio-temporal analysis in Spain and put in evidence GDP per capita and presence of mass transit systems affect the incidence of COVID-19. Paul et al. (2020) used a Bayesian hierarchical spatiotemporal model based on the Markov Chain Monte Carlo algorithm to identify hotspots of the virus diffusion. Ghosh and Cartone (2020) used several spatial indicators over different time periods to assess the spatial effects and spatiotemporal patterns of the outbreak of COVID-19 in different regions of Italy.

The goal of this research is to investigate the spatio-temporal spread of COVID-19 in Italy, by using and generalizing the main work hypothesis of Sebastiani (2020). In this paper, the author identified the provinces most affected by SARS-CoV-2 by comparing their total number of infections over a given threshold. If this number is higher than 1000, then the province is classified as critical. On April 5, 2020, more than 1000 confirmed COVID-19 cases were detected in 33 provinces; Fig. 2 shows the number of confirmed cases and the ratio of confirmed cases/resident population (incidence) of these provinces; population census data are available at https://dati.istat.it/Index.aspx?QueryId = 18460.

Fig. 2

33 Italian provinces with more than 1000 confirmed COVID-19 cases in Italy (data at 5th April 2020).

As displayed in Fig. 3 , the provinces listed in Fig. 2 are in the proximity of 4 main highways belonging to the International E-road network (United Nations Economic Commission for Europe, 2016):

Fig. 3

Italian provinces with more than 1000 positive for COVID-19 on April 5, 2020.

E35 (European route from Amsterdam, Netherlands, to Rome, Italy): its branch between Milano and Rome is labeled by the green line in Fig. 3;

E70 (European route from A Coruña, Spain, to Poti, Georgia): its branch between Torino and Venezia is given by the blue line in Fig. 3;

E55 (European route from Helsingborg, Sweden, to Kalamáta, Greece): its branch between Rimini and Ancona is the yellow line in Fig. 3;

E45 (European route from Alta, Norway, to Gela, Italy): its branch between Bolzano and Naples is the purple line in Fig. 3.

In particular, the provinces with the highest incidence of infection as of April 5, 2020 (Piacenza, Lodi and Cremona) are close or distant less than 40 km from the intersection of E35 and E70. Other relevant highway connections in the Northern Italy cannot be overlooked, namely, A7 between Milano and Genova; A6 between Torino and Cuneo; A4 between Bergamo and Brescia; A27 between Venezia and Treviso.

Here, a method based on the strategy suggested by Sebastiani (2020) is exploited to study the time evolution of the spatial distribution of COVID-19 in the Northern Italy during the first-wave scenario (Gibertoni et al., 2021). Fixed a day, rather than a single-valued threshold, a sequence of increasing thresholds is compared with the cumulative number of confirmed infections of all the provinces in the Northern Italy, yielding a classification of provinces in view of the smallest threshold overtaken by this comparison. The evolution along time of this classification is then compared with the highway network (Yao et al., 2021), to identify similarities hinting a role in the rapid spreading COVID-19 (Fig. 4 ).

Fig. 4

Number of confirmed COVID-19 cases in the northern Italian provinces at 26th February (Ministero della Salute).

On the one hand, the information concerning the transport network can be read as a set of deterministic data while, on the other, data concerning the number of infections can be viewed as observations of random variables and thus affected by error. At the beginning of the first wave, indeed, the measurement of the new infections has been affected by several kinds of errors such as, for instance, the one related to the inadequacy of the national health system to develop an efficient and fast strategy to detect new cases and to track new clusters.

It is crucial to remark that this analysis has a purely theoretical statistical character and while it aims to treat the problem by a statistical-engineering perspective, it does not aspire to be exhaustive from the medical point of view. The model focuses on the Italian epidemiological during the first wave, but it can be implemented to other time periods and areas in order to identify similarities and contrasts. However, this study does not deepen how cities can manage urban planning, mobility and land-use in the “new normal” (Corazza et al., 2021).

Methods

As aforementioned, the aim of this paper focuses on the spatio-temporal spread of COVID-19 in Italy. This work can be viewed as a generalization of the hypothesis of Sebastiani (2020), where the provinces most affected by SARS-CoV-2 were identified by comparing the cumulative number of infections over a given threshold.

By the point of view of the highway network and geographic proximity, this study has considered the following geographic characteristics:

1. neighborhood boundaries between examined territories;

2. direct highway connection.

The geographic adjacency between territories has been treated by the construction of a square adjacency matrix A. From now on, will denote the generic element of an adjacency matrix.

As far as the highway network is concerned, the presence of a highway exit within 25 km of a provincial capital has been considered as the crucial criterium to accept the hypothesis the virus could have reached the city by means of highways. This assumption has been exploited to produce a binary connection matrix C, whose entry is 1 where there exists a direct connection, otherwise it is 0. Clearly, both the matrices A and C are symmetric, that is, = and .

Finally, the third squared and symmetric matrix to be here defined is B, whose entry collects the length in km of the highway branch between the provincial capitals q and k.

The distance matrix D is then defined as the product between the three aforementioned matrices, i.e.,

Since A, B, and C are symmetric, each entry according to Equation 1:

and can be viewed as the kilometric distance between the two provinces q and k, after a thresholding procedure involving proximity and existence of a direct highway connection.

As far as the number of infections is considered, in order to guarantee a consistent mathematical approach, four assumptions on the behavior of the virus have been stated:

the virus has a “reasonable” but not “excessively high” transmissibility;

the percentages of symptomatic and “asymptomatic” individuals (i.e., infected with the virus who do not reveal any characteristic symptoms of COVID-19, but play a role in the propagation of the pandemic) with respect to the total number of patients do not depend on the geographic location of the virus, but are equally distributed over the territory;

the COVID-19 curve levels off 15 days after the start of intervention measures (i.e., lockdown efforts);

the examined data are only marginally affected by the management of the COVID-19 pandemic by different regions and local authorities.

To clarify the meaning and the role of Assumption 1, a particularly high transmission rate would not justify a spread of the virus through single outbreaks along the main transport lines, as observed in the early phase in Italy. On the other hand, this rate has to be significant because new clusters can be generated by the presence of only one infected person. This aspect is pivotal because the mathematical model used in this study assumes that transport dynamics are relevant, since even a single infected person can cause a new outbreak, but without infecting all contacted individuals.

Having regard to the time trend of the COVID-19,

1. the starting day of the analysis was set as the February 27. The time observation was set according to the space–time progression of the measures for the containment and management of the epidemiological emergency by COVID-19 adopted in Italy (Collivignarelli et al., 2020) from February 23 (11 municipalities in the Northern Italy were identified as the center of the two main Italian clusters and placed under quarantine) (Boccia et al., 2020, Calabrese et al., 2021) to March 10 (lockdown was extended to the whole territory) (Mishra et al., 2020);

Increasing varying time intervals were used to define a time sequence , where j is the number of days since the starting day: corresponds to February 27.

During the first week (), the hard-hit provinces were few, and it was easier to follow day by day the transmission and to provide then a reasonable interpretation of the data. Therefore, a daily classification of the cumulative number of new infections has been developed in order to perform a comparison with the distance matrix representing the highway connections. Moreover, in the following weeks the Italian doubling time of the COVID-19 epidemic (i.e., the number of days required for the number of cases to double, based on the rate of cumulative increase in) has been considered (La Maestra et al., 2020) to carry out this study;

3. A threshold sequence , has been defined, that is, a sequence of increasing natural numbers chosen to be iteratively compared to the cumulative number of infected in each considered province and at each value. The arbitrary threshold sequence adopted in this study is:

The maximum value (i.e., 2000) was assumed as the condition when epidemic gets out of control in the examined province.

4. Official freely available data about the number of infected were obtained from the website of the Italian Ministry of Health, where the number of infected people, day by day, split by regions and provinces (Ministero della Salute) is available; from now on the categorical index will denote the reference province, so that the natural-valued element denotes the cumulative number of infections in the province at the time .

For any province and time , and for each value of the index , comparing the number of infections and the threshold value yields a binary-response model (Horowitz and Savin, 2001) whose binary entry is given by Equation (2),

By the heuristic point of view, the indexes and identify province and day, which can be read as the “space” and “time” coordinates of the response model respectively. Thus, the sequence , , is a binary non increasing sequence, whose sum (see Equation 3) counts the number of exceeded thresholds,

Then can take values from 0 to 6; the higher is its value, the more critical is the spread of the virus in the province at the day , in order to have a synthetic outline of the SARS-CoV-2 transmission. It is thus natural to use this value to establish a classification strategy, collecting all the possible in the infection matrix whose rows sort the different provinces and columns are labeled by the different times

Fixed from the corresponding column of only the non-zero elements are selected, to build the reduced infection vector (of varying size), aimed to identify which provinces are significantly affected by the virus and “how critical” is its diffusion, since the numerical values of its entries provide a qualitative hint (the order of magnitude) on the spread of the virus.

For each chosen day , , each vector has been used to draw an undirected graph using the software R-studio(Csardi and Nepusz 2006), whose circular node sizes are proportional to the intensity

Results

In order to provide the reader with a clearer and shorter depiction of the matrix A described in the previous section, it has been split into blocks of two types: the first considers neighboring provinces within a given region (AP), while the second relates to neighboring provinces between neighboring regions (AN). As examples, Table 1, Table 2, Table 3 list the matrices AR, AP, and AN related to the Italian Northern regions (ARN), provinces of Lombardia (APL), and provinces of Lombardia and Emilia Romagna (APL,ER), respectively.

Table 1

Adjacency matrix of the Italian Northern regions (ARN).

Italian Northern region	Emilia Romagna	Friuli Venezia Giulia	Liguria	Lombardia	Piemonte	Trentino Alto Adige	Valle d’Aosta	Veneto
Emilia Romagna	0	0	1	1	1	0	0	1
Friuli Venezia Giulia	0	0	0	0	0	0	0	1
Liguria	1	0	0	0	1	0	0	0
Lombardia	1	0	0	0	1	1	0	1
Piemonte	1	0	1	1	0	0	1	0
Trentino Alto Adige	0	0	0	1	0	0	0	1
Valle d’Aosta	0	0	0	0	1	0	0	0
Veneto	1	1	0	1	0	1	0	0

Table 2

Adjacency matrix of the Lombardia provinces (APL).

Province	Bergamo	Brescia	Como	Cremona	Lecco	Lodi	Mantova	Milano	Monza e Brianza	Pavia	Sondrio	Varese
Bergamo	0	1	0	0	1	0	0	1	0	0	0	0
Brescia	1	0	0	1	0	0	0	1	0	0	0	0
Como	0	0	0	0	0	0	0	1	0	0	0	0
Cremona	0	1	0	0	0	1	1	1	0	0	0	0
Lecco	1	0	0	0	0	0	0	0	0	0	0	0
Lodi	0	0	0	1	0	0	1	1	0	0	0	0
Mantova	0	0	0	1	0	1	0	1	0	0	0	0
Milano	1	1	1	1	0	1	1	0	1	1	1	1
Monza e Brianza	0	0	0	0	0	0	0	1	0	0	1	0
Pavia	0	0	0	0	0	0	0	1	0	0	0	0
Sondrio	0	0	0	0	0	0	0	1	1	0	0	0
Varese	0	0	0	0	0	0	0	1	0	0	0	0

Table 3

Adjacency matrix of the Lombardia and Emilia Romagna (APL,ER).

		Emilia Romagna
	Province	Bologna	Ferrara	Forlì-Cesena	Modena	Parma	Piacenza	Ravenna	Reggio Emilia	Rimini
Lombardia	Bergamo	0	0	0	0	0	0	0	0	0
	Brescia	0	0	0	0	0	0	0	0	0
	Como	0	0	0	0	0	0	0	0	0
	Cremona	0	0	0	0	1	1	0	0	0
	Lecco	0	0	0	0	0	0	0	0	0
	Lodi	0	0	0	0	0	1	0	0	0
	Mantova	0	1	0	1	1	0	0	1	0
	Milano	0	0	0	0	0	0	0	0	0
	Monza-Brianza	0	0	0	0	0	0	0	0	0
	Pavia	0	0	0	0	0	1	0	0	0
	Sondrio	0	0	0	0	0	0	0	0	0
	Varese	0	0	0	0	0	0	0	0	0

For the sake of the simplicity, also in this case C will be split into two blocks, the first considering provinces of a given region (CP), and the second related to provinces of neighboring regions (CN). Table 4, Table 5 present the binary connection matrices defined for the provinces of Lombardia (CPL) and for the provinces of Lombardia and Emilia Romagna (CNL,ER), respectively.

Table 4

Connection matrix of the Lombardia provinces (CPL).

Province	Bergamo	Brescia	Como	Cremona	Lecco	Lodi	Mantova	Milano	Monza e Brianza	Pavia	Sondrio	Varese
Bergamo	0	1	1	1	1	1	1	1	1	1	0	1
Brescia	1	0	1	1	1	1	1	1	1	1	0	1
Como	1	1	0	1	0	1	1	1	0	1	0	0
Cremona	1	1	1	0	0	1	0	1	1	1	0	1
Lecco	1	1	0	0	0	0	0	0	0	0	0	0
Lodi	1	1	1	1	0	0	0	1	1	1	0	1
Mantova	1	1	1	0	0	0	0	1	1	0	0	1
Milano	1	1	1	1	0	1	1	0	1	1	0	1
Monza e Brianza	1	1	0	1	0	1	1	1	0	1	0	1
Pavia	1	1	1	1	0	1	0	1	1	0	0	1
Sondrio	0	0	0	0	0	0	0	0	0	0	0	0
Varese	1	1	0	1	0	1	1	1	1	1	0	0

Table 5

Connection matrix of the Lombardia and Emilia Romagna (CPL,ER).

		Lombardia
	Province	Bergamo	Brescia	Como	Cremona	Lecco	Lodi	Mantova	Milano	Monza e Brianza	Pavia	Sondrio	Varese
Emilia Romagna	Bologna	1	1	1	0	0	1	1	1	1	1	0	1
	Ferrara	1	1	1	1	0	1	0	1	1	1	0	1
	Forlì-Cesena	1	1	1	1	0	1	1	1	1	1	0	1
	Modena	1	1	1	1	0	1	1	1	1	1	0	1
	Parma	1	1	1	1	0	1	0	1	1	1	0	1
	Piacenza	1	1	1	1	0	1	0	1	1	1	0	1
	Ravenna	1	1	1	1	0	1	1	1	1	1	0	1
	Reggio nell’Emilia	1	1	1	1	0	1	1	1	1	1	0	1
	Rimini	1	1	1	1	0	1	1	1	1	1	0	1

Table 6 presents the distance matrix obtained for the provinces of Lombardia (DL) (values in km).

Table 6

Distance matrix of the Lombardia provinces (DL).

Province	Bergamo	Brescia	Como	Cremona	Lecco	Lodi	Mantova	Milano	Monza e Brianza	Pavia	Sondrio	Varese
Bergamo	0	53	87	103	49	65	143	60	43	90	115	99
Brescia	53	0	133	58	88	86	98	105	89	117	143	145
Como	87	133	0	148	32	91	226	50	48	88	111	55
Cremona	103	58	148	0	147	70	70	98	108	85	226	155
Lecco	49	88	32	147	0	90	176	64	42	108	80	80
Lodi	65	86	91	70	90	0	135	41	51	48	170	99
Mantova	143	98	226	70	176	135	0	195	178	151	257	235
Milano	60	105	50	98	64	41	195	0	27	45	157	58
Monza e Brianza	43	89	48	108	42	51	178	27	0	69	121	68
Pavia	90	117	88	85	108	48	151	45	69	0	188	97
Sondrio	115	143	111	226	80	170	257	157	121	188	0	159
Varese	99	145	55	155	80	99	235	58	68	97	159	0

Table 7, Table 8 present the results from thresholds implementation on infected in the Northern provinces during the surveyed period (February 27-March 28): the results refer to the first examined threshold (i.e., 50 infected in each examined province). Therefore each binary element t in Table 7, Table 8 describes the epidemic condition in the p-th province (row) at the j-th day (column).

Table 7

Threshold matrix of the Italian Northern provinces (T50).

	Sj	S1	S2	S3	S4	S5	S6	S7	S8	S9	S10	S11
	Day	27-feb	28-feb	29-feb	01-mar	02-mar	03-mar	04-mar	07-mar	12-mar	17-mar	25-mar
Region	Province
Piemonte	Alessandria	0	0	0	0	0	0	0	0	1	1	1
	Asti	0	0	0	0	0	0	0	0	1	1	1
	Biella	0	0	0	0	0	0	0	0	0	1	1
	Cuneo	0	0	0	0	0	0	0	0	0	1	1
	Novara	0	0	0	0	0	0	0	0	0	1	1
	Torino	0	0	0	0	0	0	0	1	1	1	1
	Verbania	0	0	0	0	0	0	0	0	0	1	1
	Vercelli	0	0	0	0	0	0	0	0	0	1	1
Valle d'Aosta	Aosta	0	0	0	0	0	0	0	0	0	1	1
Lombardia	Bergamo	1	1	1	1	1	1	1	1	1	1	1
	Brescia	0	0	0	0	1	1	1	1	1	1	1
	Como	0	0	0	0	0	0	0	0	1	1	1
	Cremona	1	1	1	1	1	1	1	1	1	1	1
	Lecco	0	0	0	0	0	0	0	0	1	1	1
	Lodi	1	1	1	1	1	1	1	1	1	1	1
	Mantova	0	0	0	0	0	0	0	0	1	1	1
	Milano	0	0	0	0	1	1	1	1	1	1	1
	Monza e Brianza	0	0	0	0	0	0	0	1	1	1	1
	Pavia	0	0	1	1	1	1	1	1	1	1	1
	Sondrio	0	0	0	0	0	0	0	0	0	1	1
	Varese	0	0	0	0	0	0	0	0	1	1	1
Liguria	Genova	0	0	0	0	0	0	0	0	1	1	1
	Imperia	0	0	0	0	0	0	0	0	0	1	1
	La Spezia	0	0	0	0	0	0	0	0	0	1	1
	Savona	0	0	0	0	0	0	0	0	1	1	1
Trentino Alto Adige	Bolzano	0	0	0	0	0	0	0	0	1	1	1
	Trento	0	0	0	0	0	0	0	0	1	1	1
Veneto	Belluno	0	0	0	0	0	0	0	0	0	1	1
	Padova	1	1	1	1	1	1	1	1	1	1	1
	Rovigo	0	0	0	0	0	0	0	0	0	0	1
	Treviso	0	0	0	0	1	1	1	1	1	1	1
	Venezia	0	0	0	0	0	0	1	1	1	1	1
	Verona	0	0	0	0	0	0	0	1	1	1	1
	Vicenza	0	0	0	0	0	0	0	0	1	1	1
Friuli Venezia Giulia	Gorizia	0	0	0	0	0	0	0	0	0	0	1
	Pordenone	0	0	0	0	0	0	0	0	0	1	1
	Trieste	0	0	0	0	0	0	0	0	1	1	1
	Udine	0	0	0	0	0	0	0	0	0	1	1
Emilia Romagna	Bologna	0	0	0	0	0	0	0	0	1	1	1
	Ferrara	0	0	0	0	0	0	0	0	0	1	1
	Forlì-Cesena	0	0	0	0	0	0	0	0	0	1	1
	Modena	0	0	0	0	0	0	0	1	1	1	1
	Parma	0	0	0	1	1	1	1	1	1	1	1
	Piacenza	1	1	1	1	1	1	1	1	1	1	1
	Ravenna	0	0	0	0	0	0	0	0	0	1	1
	Reggio nell’Emilia	0	0	0	0	0	0	0	0	1	1	1
	Rimini	0	0	0	0	0	0	0	1	1	1	1

Table 8

Infection matrix of the Italian Northern provinces (I).

	Sj	S1	S2	S3	S4	S5	S6	S7	S8	S9	S10	S11
	Day	27-feb	28-feb	29-feb	01-mar	02-mar	03-mar	04-mar	07-mar	12-mar	17-mar	25-mar
Region	Province
Piemonte	Alessandria	0	0	0	0	0	0	0	0	2	3	4
	Asti	0	0	0	0	0	0	0	0	1	1	3
	Biella	0	0	0	0	0	0	0	0	0	1	3
	Cuneo	0	0	0	0	0	0	0	0	0	2	3
	Novara	0	0	0	0	0	0	0	0	0	2	4
	Torino	0	0	0	0	0	0	0	1	2	4	6
	Verbania	0	0	0	0	0	0	0	0	0	1	3
	Vercelli	0	0	0	0	0	0	0	0	0	2	3
Valle d'Aosta	Aosta	0	0	0	0	0	0	0	0	0	2	3
Lombardia	Bergamo	1	2	2	3	3	3	3	4	6	6	6
	Brescia	0	0	0	0	1	1	2	3	5	6	6
	Como	0	0	0	0	0	0	0	0	1	3	4
	Cremona	1	2	2	3	3	3	3	4	5	6	6
	Lecco	0	0	0	0	0	0	0	0	2	3	5
	Lodi	2	2	3	3	3	3	4	4	5	5	5
	Mantova	0	0	0	0	0	0	0	0	2	3	5
	Milano	0	0	0	0	1	1	2	3	5	6	6
	Monza e Brianza	0	0	0	0	0	0	0	1	2	3	5
	Pavia	0	0	1	1	1	2	2	3	3	4	5
	Sondrio	0	0	0	0	0	0	0	0	0	1	3
	Varese	0	0	0	0	0	0	0	0	1	3	3
Liguria	Genova	0	0	0	0	0	0	0	0	1	3	4
	Imperia	0	0	0	0	0	0	0	0	0	1	2
	La Spezia	0	0	0	0	0	0	0	0	0	1	2
	Savona	0	0	0	0	0	0	0	0	1	1	3
Trentino Alto Adige	Bolzano	0	0	0	0	0	0	0	0	2	3	4
	Trento	0	0	0	0	0	0	0	0	2	3	5
Veneto	Belluno	0	0	0	0	0	0	0	0	0	2	3
	Padova	1	1	1	2	2	2	2	3	3	4	5
	Rovigo	0	0	0	0	0	0	0	0	0	0	1
	Treviso	0	0	0	0	1	1	1	2	3	4	5
	Venezia	0	0	0	0	0	0	1	1	3	3	4
	Verona	0	0	0	0	0	0	0	1	2	3	5
	Vicenza	0	0	0	0	0	0	0	0	2	3	4
Friuli Venezia Giulia	Gorizia	0	0	0	0	0	0	0	0	0	0	1
	Pordenone	0	0	0	0	0	0	0	0	0	1	3
	Trieste	0	0	0	0	0	0	0	0	1	2	3
	Udine	0	0	0	0	0	0	0	0	0	2	3
Emilia Romagna	Bologna	0	0	0	0	0	0	0	0	2	3	5
	Ferrara	0	0	0	0	0	0	0	0	0	1	3
	Forlì-Cesena	0	0	0	0	0	0	0	0	0	2	3
	Modena	0	0	0	0	0	0	0	1	2	3	5
	Parma	0	0	0	1	1	1	2	3	3	4	5
	Piacenza	1	1	2	2	3	3	3	3	4	5	6
	Ravenna	0	0	0	0	0	0	0	0	0	2	3
	Reggio nell’Emilia	0	0	0	0	0	0	0	0	2	3	5
	Rimini	0	0	0	0	0	0	0	2	3	4	5

The threshold matrices T100, T200, T500, T1000 and T2000 are not herein disclosed for the sake of brevity. However, they allowed compilation of the infection matrix I (Table 8).

Fig. 5a to g present the graphs obtained from qualitative infection vectors having regard to the fixed thresholds:

Fig. 5

Undirected graphs of qualitative infection vectors. a) February 27; b) February 28; c) February 29; d) March 1; e) March 2; f) March 3; g) March 4.

at day 1 (February 27) the central node, which has the largest area, is Lodi (LO), that is the province of the first Italian outbreak and it is over the first two thresholds; other provinces over-thresholds are Bergamo (BG), Cremona (CR), Piacenza (PC) (near Lodi) and Padova (PD), that is the province of the second Italian outbreak (Fig. 5a);

on the second day, the involved provinces are unchanged, but the size of the nodes area highlight the infection progression (Fig. 5b);

on the third day, a new province (Pavia-PV) is over-thresholds, while the pre-existing ones have a growing number of infected (Fig. 5c);

on day 4 (March 1), Parma (PR) joins the graphed provinces (Fig. 5d);

the next day (March 2) three provinces (Milano-MI, Brescia-BR, and Treviso-TR) reach the first threshold (i.e. 50) (Fig. 5e);

on day 6 (Fig. 5f) the number of infections increased, while the provinces remained unchanged;

on day 7 (March 4, Fig. 5g) a new province (Venezia) is in the graph.

It should be noted that over time the new provinces in Fig. 5, characterized by smaller vertices, are connected with at least one province characterized by higher vertices through direct highway connections.

Fig. 6a to d presents the graphs for days 10 (March 7), 15 (March 12), 20 (March 17), and 28 (March 25), respectively:

Fig. 6

Undirected graphs of qualitative infection vectors. a) March 7; b) March 12; c) March 17; d) March 25.

at day 10 (March 7) the most affected provinces are Lodi, Cremona, and Bergamo, adjacent to the first outbreak (Fig. 6a), while the Veneto provinces are less affected, except Padova (province of the second outbreak in Vò Euganeo);

as the days go by, both the number of over-threshold provinces and the size of the nodes increase (Fig. 6b–d): the most affected provinces tend to become centralized in the graph, while the new ones are peripheral with a tendency to centralize over time. At day 28 (Fig. 6d), Torino (TO), Milano (MI), Bergamo (BG), Brescia (BR), Cremona (CR), and Piacenza (PC) are over the upper threshold (i.e., 2000 confirmed cases of COVID-19), while the province of Rovigo (RO) appears for the first time in the graphs of qualitative infection. Therefore, at March 25 all the provinces in the northern Italy (except to Gorizia - GO) are above at least one threshold.

Discussion

This study investigated the influence of key built environment factors such as interurban mobility and social contacts on the spread of the novel coronavirus disease in Italy during the first wave (Galiano et al., 2021). However, the proposed model can be applied to any other study areas to identify the relationship between the high-speed road network and the disease transmission in different contexts and time periods (Cheshmehzangi et al., 2021). Indeed, in the pre-COVID-19 the literature demonstrated different transport modes (Findlater and Bogoch, 2018) as a direct factor affecting infectious diseases transmission (Zhang et al., 2020a, Zhang et al., 2020b, Peak et al., 2018, Arthur et al., 2017, López-Quílez, 2019). The recent COVID-19 literature focused on built environment factors (Ma et al., 2021), urban attributes (Barak et al., 2021), and living conditions (Sahasranaman and Jensen, 2021) to investigate the income on cases growth. The built environment plays a pivotal role in people health and well-being (Angiello, 2021): density of schools, commercial and sport centers, proportion of built-up areas, building systems and nighttime light, socioeconomic and demographic characteristics have been identified have as key drivers of the SARS-CoV-2 transmission and infection (Verma et al., 2021, Barak et al., 2021, Corazza et al., 2021). Some studies proposed temporal and spatial analysis of COVID-19 spreading focusing on air flows and dispersion of aerosolized virus in polluted areas (Zheng et al., 2021), meteorological variables (Wang et al., 2021), traffic volumes (Mu et al., 2020). Network analyses have been carried out on cross-country pandemic air connectedness (Chu et al., 2021), maritime transport (Wang et al., 2022), high-speed railway (Zhang et al., 2020a, Zhang et al., 2020b): the attention was paid on mass public transport. Conversely, in this study the authors investigated how the highway transportation network in the Northern Italy exacerbated the large-scale spreading of COVID-19 during the first wave. This study gives the private transport a pivotal role in the Italian scenario: such conclusion requires further analyses in terms of traffic flow and comparison with other States to confirm the results and identify critical issues. Time-series plots and spatial–temporal maps allowed the analysis of the outbreak in 8 regions and 47 provinces and the results confirmed that mobility was a significant direct factor to the COVID-19 spread. Indeed, the comparison between the graphs in Fig. 5, Fig. 6 and the map of the highway network highlight a correspondence both with the binary adjacency matrices and the connection matrices:

having regard to E35, the virus moved along the highway starting from the Lombardia outbreak of the Lodi province North to Milano and South to Bologna (Fig. 7 ). In the first four days the spread of the virus follows the highway from Lodi southwards to Parma via Piacenza. Then the Milano province is added to day 5, the Modena province is added to day 10, finally Bologna and Reggio nell’Emilia provinces are added to day 13. The spread southwards has a delay in succession, because Reggio nell’Emilia reaches the first threshold (i.e. 50) a couple of days after Modena: the geographical distances (34 km via E35, 27 km via SS9) and the compact and continuous urbanization make interchangeable these two cities, therefore the observed delay does not have significance from a statistical point of view. Moreover, it should be noted that the daily analysis of the qualitative infection vector revealed that the Bologna province had a faster epidemic progression than the Reggio nell’Emilia province, probably due to railway hub in Bologna. In the analyzed territory the railway network overlaps E35 and allows a greater number of connections between the provinces: the railway system may have bypassed the highway network and the intermediate provinces along the highway.

Fig. 7

Virus spread along E35.

having regard to E70, from Piacenza and Cremona (i.e. provinces affected by confirmed COVID-19 cases at day 1) the virus spread both to the East and to the West. It reaches Brescia at day 5, Torino at day 10, while Alessandria and Asti are added to day 13. In this case, as for Bologna along E35, it is necessary to suppose a pivotal role of the high speed railway line that between Torino and Milano (See Fig. 8 ).

Fig. 8

Virus spread along E70.

on the other hand, E70 played a role in the virus spread from the second Italian outbreak (in Vò Euganeo Padova province, in the Veneto region) to the Friuli Venezia Giulia region. Indeed, from Padova it headed towards Trieste and Udine passing through Treviso and Venezia along E70 (Fig. 9 ). In fact, on day 1 only the Padova province is over-threshold (Fig. 5a), and after day 7 other two Venetian provinces are added (Fig. 5g): Venezia and Treviso; on day 15 the Trieste province joins to the qualitative infection graph (Fig. 6).

Fig. 9

Virus spread along E70.

having regard to E45, Trento and Bolzano provinces are over-threshold since day 15 (Fig. 10 ), whereas Verona province was over-threshold since day 10 (Table 8);

Fig. 10

Virus spread along E45.

having regard to E55 (Fig. 11 ), the spatial–temporal spread of COVID-19 confirms the explored hypothesis: on day 20 (i.e., after the Venezia, Treviso, and Trieste provinces are over-threshold) the Udine province joins to the qualitative infection graph (Fig. 6). Therefore, the layout of the highway network could justify the Udine’s delay compared to Trieste. Moreover, it could be justified by the transportation role of Trieste, whose role as European hub for maritime transport may have speed up the spread.

Fig. 11

Virus spread along E55.

The discussed examples could not demonstrate the overall proposed model, but they are samples to evaluate how the model works well in reconstructing the spread of COVID-19 infection along the Italian highway transport network. Indeed, a province that needs a separate discussion is that of Bergamo. It contains the first large municipality (Bergamo) that is near Lodi and at already day 1 it had a high number of infected. The trend observed in the city of Bergamo can be interpreted as an exception: its province became the starting point of secondary outbreaks, or the third outbreak as Lodi and Padova. Indeed, Bergamo acted as a link between the Lombard (i.e. Lodi province) and the Venetian (i.e. Padova province) outbreak: the virus arrived in Verona both from Bergamo passing through Brescia and through the Paduan outbreak, along both directions of E70 (Fig. 12 ): on day 1 both the Bergamo and Padova provinces are over-threshold, the Brescia province joins at day 5, and Verona at day 10.

Fig. 12

Virus spread along E55.

Conclusions

Starting from previous analyses in the literature that believe the virus traveled on roads, this paper proposes a model to investigate the spatio-temporal pattern of COVID-19 in Italy during the first-wave in Spring 2020. Since 24th February to 28th March the number of people confirmed having contracted the coronavirus SARS-CoV-2 in the northern provinces has been monitored with a variable time interval. The analysis was conducted day-by-day in the first week, when few provinces were involved, in order to conduct a thorough investigation. Given the large number of involved provinces after the first week, the time observation was set according to the doubling time of the infected in order to summarize the results (3 days between March 1 and 10, 5 days between March 11 and 20, 10 days in the rest). The doubling concept was implemented to six thresholds of infected to be compared to the number of confirmed cases in each province: thresholds approximate a geometric progression (i.e. 50, 100, 200, 500, 1000, and 2000 confirmed infected people) to have a detailed survey of the infection process. The comparison between the number of infected and the assumed thresholds allowed definition of binary matrices and vectors to describe the temporal dynamics of the virus. Finally, a binary model of the highways belonging to the International E-road network was defined to interpret the spatial spread of the COVID-19 epidemic in Italy. The well-fitting model results demonstrated that the highway network contributed to the diffusion of COVID-19 during the first wave in the northern Italy. From the first two outbreaks, the virus regularly moved along E35 to the Milano and Bologna provinces, along E70 to the Torino and Trieste provinces, along E45 to the Trento and Bolzano provinces, and along E55 to the Udine province. Moreover, the city of Bergamo acted as a link between the Lombard and the Venetian outbreaks. Observed delays in succession can be caused by different transport mode (e.g., railway or sea) or characteristics of the built environment that may have speed up the spread.

This approach can be implemented to study more recent waves within Italy under the hypothesis that private mobility along highways had the same key role in different time periods. The method is useful to interpret the epidemiological data collected in other countries: the spread of the infections over countries could be connected by several means of transport (both private and public). In particular, the Italian and non-Italian studies could help to identify variables other than high-speed road network (e.g., built environment, socio-economical fabric, public health recommendations) that affected the COVID-19 diffusion over time and space. Furthermore, in our opinion, this kind of approach could be easily applied as a preliminary study to any diffusion process. Indeed, it aims to provide a qualitative study of its growth rate and in this sense we consider it very versatile.

CRediT authorship contribution statement

Marco De Angelis: Data curation, Investigation. Claudio Durastanti: Conceptualization, Methodology, Writing – original draft, Writing – review & editing, Visualization. Matteo Giovannoni: Data curation, Investigation. Laura Moretti: Conceptualization, Methodology, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Region	Province
Piemonte	Alessandria
	Asti
	Biella
	Cuneo
	Novara
	Torino
	Verbania
	Vercelli
Valle d'Aosta	Aosta
Lombardia	Bergamo
	Brescia
	Como
	Cremona
	Lecco
	Lodi
	Mantova
	Milano
	Monza e Brianza
	Pavia
	Sondrio
	Varese
Liguria	Genova
	Imperia
	La Spezia
	Savona
Trentino Alto Adige	Bolzano
	Trento
Veneto	Belluno
	Padova
	Rovigo
	Treviso
	Venezia
	Verona
	Vicenza
Friuli Venezia Giulia	Gorizia
	Pordenone
	Trieste
	Udine
Emilia Romagna	Bologna
	Ferrara
	Forlì-Cesena
	Modena
	Parma
	Piacenza
	Ravenna
	Reggio Emilia
	Rimini