#### Abstract

During the COVID-19 epidemic, draconian countermeasures forbidding nonessential human activities have been adopted in several countries worldwide, providing an unprecedented setup for testing and quantifying the current impact of humankind on climate and for driving potential sustainability policies in the postpandemic era from a perspective of complex systems. In this study, we consider heterogeneous sources of environmental and human activity observables, considered as components of a complex socioenvironmental system, and apply information theory, network science, and Bayesian inference to analyze their structural relations and nonlinear dynamics between January 2019 and August 2020 in northern Italy, i.e., before, during, and after the national lockdown. The topological structure of a complex system strongly impacts its collective behavior; therefore, mapping this structure is essential to fully understand the functions of the system as a whole and its fragility to unexpected disruptions or shocks. To this aim, we unravel the causal relationships between the 16 environmental conditions and human activity variables, mapping the backbone of the complex interplay between intervening physical observables—such as NO_{2} emissions, energy consumption, intervening climate variables, and different flavors of human mobility flows—to a causal network model. To identify a tipping point during the period of observation, denoting the presence of a regime shift between distinct network states (i.e., before and during the shock), we introduce a novel information-theoretic method based on statistical divergence widely used in statistical physics. We find that despite a measurable decrease in NO_{2} concentration, due to an overall decrease in human activities, locking down a region as a climate change mitigation is an insufficient remedy to reduce emissions. Our results provide a functional characterization of socioenvironmental interdependent systems, and our analytical framework can be used, more generally, to characterize environmental changes and their interdependencies using statistical physics.

#### 1. Introduction

Panels show the time series for each of the 16 variables considered in this study; blue and red curves correspond to the time courses for years 2019 and 2020, respectively. The light-red band corresponds to the lockdown period in 2020, while the light-blue band corresponds to the same period in 2019. The public mobility data are available only for the year 2020.

Complex systems consist of interconnected units, which are characterized by nonlinear dynamics at the microscopic scale that lead, at larger scales, to emergent collective phenomena. Often, those units are complex subsystems interacting with each other, exhibiting a multilayer [1] or interdependent organization [2]. Complex networks provide an abstract representation for the backbone of such systems. In fact, it has been shown that shocks in one node or in one subsystem can quickly propagate to the rest of the network, driving the overall system to a catastrophic collapse [3–5]. However, understanding the resilience of a complex system can be even more challenging when its backbone cannot be directly observed and must be inferred from indirect observations, such as the time course of a set of physical observables. In fact, the collective behavior of the system is driven by the interplay between the dynamics at each component and the hidden interactions among them. In statistical physics, the reconstruction of the system’s interactions from measuring time series related to its components is usually referred to as inverse problem [6]. In the last decades, a great effort has been devoted to the solution of this problem [7–11], which has relevant implications in various disciplines, from neuroscience [12], to ecology [13], including finance [14], biology [15], and climatology [16]. This is also the case of the complex nexus between environment and human activity, where a direct relationship between environmental, climate, and human-dependent variables cannot be measured. However, mapping the functional backbone of these interdependent subsystems, what in the following we name the socioenvironmental system, is of paramount importance to understand their response, and consequently their resilience, to exogenous and endogenous perturbations.

To this aim, we focus on a specific case study, namely, the reduction of pollutants in the northern Italy, observed between March and July 2020 in response to draconian interventions due to the COVID-19 pandemic. In fact, during the first months of 2020, Italy has been one of the countries mostly affected by COVID-19 [17–25]. To prevent the spread of the SARS-CoV-2 virus, Italy adopted significant nonpharmaceutical interventions [26–29], including locking down the entire country from of March to of May. The forced closure of schools, public facilities, and places of employment drastically reduced the vehicle traffic and industrial activities, with the most relevant effects in the Lombardia region, in the northern Italy, which is also the Italian region most plagued by the virus. On the one hand, this unprecedented situation triggered a cascade of public health, social, behavioral, and economic challenges that will require years to recover. On the other hand, from a scientific perspective, one can investigate the effects of such dramatic regional and subregional interventions on the environment. In practice, the spread of COVID-19 allows one to better understand the causal and functional relations between the components of the socioenvironmental system, also shedding light on the potential effects that large-scale interventions, such as lockdowns, may have on the system’s collective behavior. This insight would also be a valuable resource to design informed policies that can be adopted to mitigate the climate crisis. To this aim, we consider this situation as a global experiment, giving us a unique opportunity to investigate the complex interactions between human activities and the environment with a particular focus on air quality, using the northern Italy as a case study where the availability of heterogeneous data sources allows one to perform a more integrative analysis. The interest in such a relationship has exploded worldwide during the lockdown period, and many recent studies focus on the reduction of atmospheric pollutants due to social restrictions in different countries [30–44]. However, some works highlight that reductions in human mobility and in industrial activity would not be sufficient to reduce air pollution, especially when meteorology is unfavorable [45]. In fact, it must be noticed that air quality does not only depend on emissions (and consequently concentrations) of pollutants but it is also strongly influenced by meteorology, which plays a crucial role in the Po valley, where the Lombardia region is located. Moreover, the effect of the topography, with the Alps to the north, contributes to making this region one of the most polluted areas around Europe.

Here, our purpose is to unravel the complex interactions between environmental conditions and human activity, taking the Lombardia region during the Italian lockdown as a case study. Mapping the network structure of a socioenvironmental system is crucial to understand its collective behavior and to acknowledge the emergence of unexpected outcomes [46, 47]. To this aim, we reconstruct the network of functional and causal relations between the observables and subsequently identify the (causal) influence that an external intervention—as the lockdown—exerts on the system. In particular, we aim at evaluating the impact of relaxing nonessential human activities, i.e., those activities not directly related to the supply of goods and commodities, on tropospheric concentrations during the lockdown. We employ representative data for environmental conditions, i.e., air pollution and meteorological conditions data, and human activity, i.e., mobility and energy load data—over a survey period of four months—February, March, April, and May—both in 2019 and in 2020, for a total of 16 variables Figure 1.

Figure 2 shows the 5-day moving average concentration of the in the two years of reference. The dashed vertical line indicates the beginning of the lockdown period on March 9, 2020.

By taking into account all these variables, we reduce the possible disturbances due to latent confounders. Afterward, we evaluated the differences in the dynamics of such variables between the period related to a situation with extremely reduced activity and baseline periods. Since the lockdown imposed social distancing and, consequently, drastically reduced human mobility, we referred to an indicator of air quality, which is sensitive to mobility changes, i.e., the nitrogen dioxide concentrations, which proves to highly depend on emissions from transportation means and industrial activity [48]. Concerning the meteorological conditions, we considered the daily average of 6 variables, i.e., surface pressure, surface net solar radiation, temperature, total precipitation, wind direction, and wind speed. We used the data made publicly available by Google [49], Apple [50], and Facebook [51] during the survey period to investigate human mobility changes. In the Lombardia region, energy consumption is mainly attributable to industry; consequently, we have considered the daily average electricity system’s total load, as a proxy of industrial activity (see Appendix A for details). From a methodological point of view, we evaluated the significance and the concomitance of variations in human activities and in environmental conditions—focusing on concentrations—in the survey period 2019–2020, using statistical and causal analysis. Through statistical tests and effect size measures [52, 53], we were able to assess the relevant variations in the time series of the considered variables. To investigate the complex nexus between the 16 variables, we leveraged on the partial correlation coefficient (PCC) [54] and Granger causality (GC) [55]. These two measures provide functional and causal topological maps of the complex system. Finally, we assess the causal impact of the lockdown intervention by taking inspiration from a relatively recent Bayesian technique, based on a state-space model, used to infer the causal impact of an advertising campaign on market sales [56]. By specifying which period in the data should be used for training the state-space model (preintervention period) and which period for computing a counterfactual prediction, this technique assesses the impact of the attributable intervention [56].

#### 2. Identifying Tipping Points in Empirical Observations

Figure 1 provides an overview of the time course of each observable used during the survey period. In Figure 2, the variation is reported in detail, comparing the two years of reference, suggesting that the average concentration during the lockdown period in 2020 is smaller than the one during the same period in 2019 (see Appendix B for the exact start and end dates of the periods used in the statistical tests). In the following, this observation will be corroborated by statistical analysis. It is worth noticing that, in general, at the end of the winter season, the concentration is expected to decrease due to the reduction in the domestic, civil, and industrial heating concurrently with the rising temperatures. In 2020, the reduction appears to be more abrupt with the beginning of the lockdown and the concentration levels remain lower for the rest of the period. The variability of the concentration during the lockdown is also lower than the one in the same period in 2019 (the Fligner test on the difference of variances [57] gives a *p* value ), and this may be due to the reduction in the fast-changing pressure variables, such as transportation. The results of *t*-tests [58] support the visual assessment of Figure 2, confirming that even though the concentrations in the prelockdown periods (in 2019 and in 2020) are statistically equivalent (the hypothesis of equal average concentration for the two periods cannot be rejected with a *t*-test *p* value ), the reduction observed during the 2020 lockdown with respect to the same period in 2019 is statistically significant (the average concentration in the two periods are different, with a *t*-test *p* value ). These results confirm that the average concentration in 2020 was significantly lower than in 2019. Now we want to evaluate also the magnitude of the observed differences by computing their effect size. Two commonly used measures of the effect size are the Cliff- and the C.L.E.S. (common language effect size) [52, 59] (see Appendix C for details). For the prelockdown periods, these measures indicate a very low effect size ( and ), whereas it is very high for the postlockdown periods ( and ). The CLES measure provides an immediate interpretation of the results, and in particular, a means that the probability that an observation from the 2020 lockdown period returns a lower value of concentration with respect to the same period in 2019 is more than 80%. A complete summary of the statistical tests is reported in Table 1.

Statistical tests on meteorological variables exclude the possibility that the average meteorological conditions in the lockdown period of the two years 2019–2020 are different (see Table 2 for detailed results). We can, therefore, hypothesize that the environmental conditions were similar and that the variation in the concentration was driven by the changes in human activities. This surmise, which considers only the average conditions, is further investigated through the causal impact analysis (see the next section).

Finally, the total load in the lockdown period in the year 2020 is significantly smaller than that recorded in the year 2019 (*t*-test on average loads difference gives a *p* value ), while there is no significant difference in the prelockdown periods.

It is not possible to test the differences in mobility trends due to the lack of publicly available data for the year 2019, but a visual analysis of Figure 1 makes evident the pronounced decrease in mobility during the lockdown.

#### 3. Building the Causal Nexus

Blue edges represent the negative partial correlation, while red edges represent positive partial correlation; the thickness of the edges is proportional to their part.

Note that only the meteorological regressors can be used for the counterfactual prediction, since they are the only variables not influenced by the lockdown intervention.

In the physics of complex systems, the reconstruction of the topological network structure of a system is a known problem attracting increasing interest in many fields [10, 11, 60–62]. Here, we apply two techniques, namely, the partial correlation and the Granger causality, to reconstruct the network of statistical and causal relations between the components of the socioenvironmental system. In Figure 3, the structures of the two resulting networks are reported, and details on partial correlation and Granger tests are shown in Figure 4 and Appendix D. The partial correlation network in Figure 3(a) shows the statistical relation between pairs of system components, after removing the effect of the other components. This structure reveals the stronger relations between all human activity variables, which synchronously vary, driven by the lockdown interventions. In particular, the variable “change in movement” is the most connected. The negative relation between “residential” and “work places” is well captured by the method, likewise for the “total precipitation,” “surface radiation,” and “temperature.” The variable is related just to the solar radiation, which may be interpreted as the influence of the seasonal effects mentioned in the introduction.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

Figure 3(b) shows the Granger causality network that represents the (Granger) causal influence of a variable onto another (see Appendix C for more detail on the definition). This network, that is directed and appears denser than the previous one, points out the influence of human activities on the concentration and also the possible influence of meteorological conditions that could cause a variation on air pollution, such as precipitation. Meaningful meteorological relations can be found in this network, whereas the human activities constitute a dense separated cluster.

On the one hand, the partial correlation and the Granger causality reveal the complex network of interaction between the variables of the system. The topology of the interaction strongly impacts the collective dynamics of the variables and is essential to understand the function of the entire system. On the other hand, for the evaluation of the impact that a large-scale “intervention”—as the lockdown restrictions—may have on the socioenvironmental system, we need to assess the causal relation between the specific event of the lockdown and the measured effects on the variables under study. In particular, we want to evaluate the causal relation between the concentration of pollutants and the lockdown intervention. To this aim, we employed a Bayesian structural time series model to build a counterfactual prediction of what would have happened to the concentration if the lockdowns were not implemented (see Appendix C for details).

The results are shown in Figures 3(c) and 3(d). The plot at the top of Figure 3(c) shows the concentration data and the counterfactual prediction for the lockdown period. The concentration during the lockdown period had an average value of mol/ , whereas in the absence of an intervention, we would have expected an average response of mol/ . The difference between concentration data and the corresponding counterfactual prediction middle plot in Figure 3(c) yields an estimate of the causal effect of the lockdown. This effect is mol/ with a interval of . In relative terms, the response variable showed a decrease of with a interval of . To obtain the cumulative impact of the lockdown, the concentration data are summed up the bottom plot in Figure 3(c) obtaining a cumulative concentration of mol/ , which would have been equal to mol/ in the absence of the lockdown. Even though the results seem to reveal a clear causal impact of the lockdown on the concentration, the computed probability of obtaining this effect by chance is . This is an evident signal of a lack of statistical significance. To sum up, considering also the results of the statistical tests described above, the reduction in the concentration is significant, but this reduction may not be caused just by curbed human activities. This fact, although restricted to the specific case study, may have relevant implications from both scientific and policy-making standpoints discussed in the following section. The color of the entries indicates the *p* value estimates. Black crosses are placed where the null hypothesis (no relation between variable) is rejected. The significant relations are depicted in the networks of Figures 3(a) and 3(b).

#### 4. Conclusions and Outlook

To understand the impact of extensive human interventions on the environment, the social and environmental systems cannot be regarded as “closed systems.” In fact, the collective behavior of the socioenvironmental system depends on the mutual feedbacks between these two counterparts, which in turn give rise to emergent properties of different nature [46]. To recognize such feedbacks, the physics of complex systems leverages on the characterization of the structure of relations between the dynamical components of a system. Capitalizing on powerful tools from network science, we revealed the structure of functional relations between the components of a socioenvironmental system, here considered as a complex network. Moreover, by coupling data analysis and causal inference we studied the impact of a large-scale intervention on such a complex system, identifying the shifts in its behavior and providing a formal identification of the related causes. To this aim, we have fused heterogeneous data sources—including human mobility, total energy load, meteorological conditions, and concentrations—in four different periods over 2019 and 2020. These variables are considered as signals coming from the different components of the socioenvironmental system. From this information, we determined the impact of the lockdown intervention on the components detecting a regime shift and assessing the significance of any different conditions by means of statistical tests. In particular, we designed and applied an information-theoretic technique to find the date of a possible regime shift in the data, which turns out to be a few days after the imposition of the lockdown. We have found evidence that, concomitantly with the reduction in both the mobility and the energy demand, the NO_{2} average concentration significantly decreases in the 2020 lockdown w.r.t the same period in the previous year. The lower variance of during the 2020 lockdown is even more visible than the decrease in the average concentration. This is attributable to the complex nature of the human-environmental system where internal factors may act as a filter on the rapid variations of the . Shedding light on this factor would unveil new possible features of the system under consideration.

The functional organization of the complex system was reconstructed from the data using partial correlation and Granger causality. The results provide a topological map of the socioenvironmental system, showing the mutual influence of the variables. Finally, we investigated the influence of the lockdown intervention on the variables relying on a counterfactual Bayesian model. Overall, the analysis detected the sign of a causal relationship between the relaxation of a broad spectrum of human activities and air pollution abatement during the lockdown in northern Italy. However, the statistical significance of the causal impact result is questionable; this has important implication from both statistical and policy-making perspectives: the lack of significance suggests that the effect might be due to the chance with non-negligible probability, but this can also be due to other uncontrollable issues related to the data, and it cannot be undoubtedly interpreted as a lack of causal relation. For example, although the meteorological-seasonal variables, which are not affected by the lockdown intervention and are related to the concentrations, are good candidates as explanatory variables of the in the Bayesian model, they may still not have sufficient predictive power on the . In fact, the ideal experimental setup should involve a stronger driver of concentration as an additional explanatory variable, which must not be affected by lockdown interventions to produce the counterfactual prediction (as further explained in Appendix C). However, it is difficult to identify such a variable, since the main sources of are anthropogenic, and virtually any anthropogenic source of the pollutant was affected by the lockdown, making these variables unusable in the causal model. From the policy-making perspective, the lack of a clear causal relation between the lockdown and the concentration reduction should point out that reducing emissions from traffic and power plants maybe not as effective as previously thought. For example, heating and cooking systems, and industrial and agricultural emissions should be better monitored and controlled. Nevertheless, our result is case specific, and its robustness should be further assessed by performing analyses on a broader geographical region. In fact, as we showed in previous sections, particular environmental conditions (e.g., meteorology and topography) have a decisive impact on pollutant concentrations. For this reason, policies should be designed assuming a systemic point of view while keeping a context-specific insight.

In conclusion, our results suggest that a lockdown intervention may reduce the pollution concentration, but it may not be a resolutive nor a definitive remedy. The weak causal impact of the intervention can be the sign of the predominant role played by the emissions from essential and unstoppable human activities. Furthermore, other external conditions (e.g., meteorology and topography) may antagonistically play, undermining the positive effects of the restrictions. In addition, living in lockdown conditions is proving to be economically unsustainable [63, 64]. In the light of these findings, we consider our approach to be indicative, but not definitive, for investigating the nexus between environmental conditions and human activity during COVID-19 lockdown in Italy. Further developments should complement our analysis with mobility data of 2019 and including more detailed data such as stratified transport (i.e., heavy and light transport with differentiated emissions). In fact, it is known that motorway traffic has considerably decreased during the lockdown, but the reduction in the circulation of heavy transport has not been so similarly comparable, due to the delivery of essential commodities. For example, in March 2020, even though the highway total vehicles dropped by , the heavy vehicles decreased by just : the latter pollutes four times more than light transport, according to the Italian National Autonomous Roads Corporation (ANAS). While the continuity in the essentials’ supply chain proves the efficiency of the region in providing indispensable services during emergencies such as the COVID-19 pandemic, it also testifies that the backbone of human activities has never really stopped. Consequently, it is plausible to hypothesize that such a backbone might play a weightier role than a nonessential one, posing serious challenges to policy- and decision-making to build a sustainable society in response to climate change.

We firmly believe that a paradigm shift toward a systemic view is necessary and fundamental, when studying the complex relations between society and the environment. This broad view has been successfully applied in this context and also on a finer scale, to understand the complex and highly nonlinear dynamics of pollutants due to their own interactions taking advantage of the lockdown conditions [65, 66]. Such an approach would be of tremendous help for decision-making processes allowing for more informed and integrated choices, especially for the development of mitigation policies in accordance with climatic and environmental goals (e.g., Sustainable Development Goals of Agenda 2030).

In this study, we have provided a practical computational framework to analyze the time course of correlated social and environmental observables, showing how network science and statistical analysis can be used to identify a regime shift and, more importantly, analyze the impact of systemic shocks such as lockdown interventions during a pandemic.

#### Appendix

#### Appendix A. Overview of the Datasets

In this work, we relied on data of nitrogen dioxide ( concentrations) from Copernicus Sentinel-5P satellite from January 1, 2019, to June 1, 2020 (TROPOMI Level 2 total column products, version 01, European Space Agency [67]). In particular, we referred to high-resolution daily concentrations of the tropospheric over the Lombardia region. Data of meteorological conditions are retrieved from the Copernicus Climate Change Service [68] (ERA5-Land reanalysis dataset) and consist of hourly data of 6 variables, i.e., surface pressure, surface net solar radiation, temperature, total precipitation, wind direction, and wind speed, over the Lombardia region. These data were averaged daily. Google mobility data [49] are provided in terms of daily length of stay at different places, e.g., residence, grocery, and parks, aggregated at the regional level. Apple data [50] are provided in terms of variation in the volume of *driving* direction requests with Facebook data [51] in terms of positive or negative *change in movement* relative to baseline (February 2020). We considered the case of the Lombardia region collectively covering the time period from January 13 to July 27, 2020. All the mobility data are available only for the year 2020. For what concerns energy load, we considered data of northern Italy for the period from January 1, 2019, to July 27, 2020, from the Italian transmission system operator *Terna* [69]. “Northern Italy” is the smallest available space aggregation, which includes Lombardia; since Lombardia has the highest energy demand in this area, northern Italy data are deemed appropriate for our purposes. The main descriptive statistics of all the variables are reported in Table 3.

#### Appendix B. Jensen–Shannon Divergence and Period Selection

To determine whether the data reflect the regime shift of the lockdown, we implemented a shifting point detection technique based on an information-theoretic measure of similarity, called Jensen-Shannon divergence. The Jensen-Shannon divergence (JSD) is an information-theoretic measure, which quantifies the difference between two probability distributions. It has the main advantage of being symmetric and that its square root defines a metric; therefore, it can be used as a distance measure to quantify (dis)similarity.

Given two discrete distributions, and Q(x), over the same probability space— in our case—one can define the Kullback-Leibler divergence (KLD) from to as

This measure is also known as relative entropy, and as by definition, it is not symmetric. It measures the amount of bits that one gains about using to model the distribution ; in fact, when then the divergence is zero bits. This relative entropy has also the disadvantage of not being upper bounded; for this reason, other measures like the JSD are also widely used.

By introducing the mixture distribution , the JSD can be defined in terms of the KLD as

It can be shown that bits. In our study, the role of and is played by two subperiods of observation; given a time series , it is split into for and for , being and the initial and final observation time, respectively. Here, represents the instance of time at which the observation is split and we apply the method explained in the text to find its value by means of the JSD.

Through this data-driven procedure, we found the date that splits the survey period into prelockdown and lockdown periods. This is used both for statistical tests and for the Bayesian state-space model. In particular, we applied the Jensen-Shannon divergence to shift subsets of the time series (of length ) for the year 2020.

The result of the process is the date of the tipping point in which the dynamics meet a regime shift. The date of the “information-theoretic” lockdown turns out to be the of March (*p* value = ), 5 days after the institutional lockdown. This result is consistent with the physical behavior of the , which has a typical lifetime of few days. Consequently, we grouped the time series in the 4 subsets used for testing:(i)Group 1.1 (prelockdown) from February 1, 2019, to March 14, 2019(ii)Group 1.2 (lockdown) from March 14, 2019, to May 5, 2019(iii)Group 2.1 (prelockdown) from February 1, 2020, to March 14, 2020(iv)Group 2.2 (lockdown) from March 14, 2020, to May 5, 2020

#### Appendix C: Statistical and Causal Analyses

To rigorously assess the differences in the time series of the considered variables, we tested the differences in means between each group (as defined in *Appendix* B) with *t*-tests and surrogate data tests. In addition, we evaluated the differences in the variances between the groups with the Fligner-Killeen test, which is a nonparametric test for homogeneity of group variances based on ranks, robust against non-normal data [57].

The magnitude of the differences observed in the time series is evaluated through two effect size measures: the Cliff- [59] and the C.L.E.S. (common language effect size) [52]. The Cliff- is computed by enumerating the number of occurrences of an observation from one group, e.g., having a higher response value than an observation from the second group (e.g., observation ), and the number of occurrences of the reverse as follows:where

The indices 1 and 2 refer to the two groups under consideration (as defined previously in this section), and the two time series and are of size and , respectively. This statistic measures the tendency of each in group 1 to be higher than each , in group 2, and it is not dependent on any assumptions whatsoever. The values of the statistic can run from (nonoverlapping distributions with smaller ) to 1.0 (nonoverlapping with smaller ). Similarly, the C.L.E.S. is defined as the probability that a randomly selected individual from one group has a higher score on a variable than a randomly selected individual from another group. This measure has the advantage of being simply interpreted as the probability of a value from one time series to be higher than a value in the other [52]. We numerically computed this measure with a brute-force approach, by repeatedly random sampling the values from the time series of the two years and computing the frequency with which the values of the first series were higher than the ones in the other.

Afterward, to obtain a statistical indication of the possible causal relations between the 16 variables, we build the partial correlation matrix. The reconstruction of network structures using correlation measures is a very well-known technique in network science (see, e.g., [10, 70, 71]). The partial correlation measures the strength and the direction of the (rank) dependence of two variables from a set of random variables when the influence of the remaining variables is removed. More precisely, we computed the partial correlation coefficient for each pair of variables and with the effects of the remaining variables removed. This can be done by performing two separate linear regressions on the variables and using the vector as a regressor:where is the vector of the regression coefficient, and is its optimal OLS estimation, after calculating the residuals and for the two variables as

The partial correlation coefficient is computed as the Pearson correlation between the residuals:

Subsequently, we tested the significance of the obtained partial correlation values comparing them to the results from a surrogate analysis, consisting in the estimation of the partial correlations using numerous reshuffled (not correlated) versions of the original time series [72]. Therefore, only the statistically significant relations were considered in the reconstructed network.

We further investigated the structural features of the system, studying the network of causal dependence of the variables through the “predictive causality” measure by C. Granger [55]. The Granger causality is a tool of increasing interest in the study of the structural features of complex systems [55, 73–76]. It is defined as a causality test between two time series and , which detects the ability of one variable to predict the other. In other words, it says that if the information about the trajectory of improves the prediction of the trajectory of , then causes . In particular, an autoregressive model is estimated for the variables and in the formwhere is the lag, and the coefficient of the model is optimized (max-likelihood). Then, if then it can be concluded that Granger causes . To obtain a statistically significant conclusion, a Wald test [77] is performed, which evaluates the null hypothesis . It is to be noted that, since in general , also the direction of the causal relation can be assigned, allowing to distinguish the situation in which Granger causes from Granger causes ; consequently, the network will be directed as well. Since the Granger test relies on an autoregressive model, the time series should be stationary. In our case, the observed time series is not stationary—according to the Augmented Dickey-Fuller test. Therefore, we differentiated the data obtaining a set of time series that passed the unit root test.

Finally, to assess the impact of the lockdown on the pollutant concentration we employed a Bayesian modeling technique for casual inference, which uses a structural time series model to predict what would have been the system evolution after an intervention, if the intervention had never occurred [56]. Structural time series models are state-space models for time series data. They can be generally defined as follows:where and are independent of all other unknowns. Equation (C7) serves as a link between the observed data and the latent state vector , which is a -dimensional vector. In equation (C8), the evolution of the latent state is described. The term is a -dimensional output vector, is a transition matrix, is a control matrix, is a scalar observation error with noise variance , and is a -dimensional system error with a state diffusion matrix where . In our case, the observation term is the vector of observation of the. The general structure of equations (C7) and (C8) can be adapted to describe different behaviors of the latent state, including local trends, seasonality, and the influence of covariates (for more details on the model structure see, e.g., [56, 78]). The covariates, which are important actors for the aim of our analysis, can be included in the model considering a static regression by setting and or by dynamical regression:where is the index for each covariate, , is the coefficient for the -th control series, and is the standard deviation of its associated random walk. This regression component can be turned into the same structure as equations (C7) and (C8) by setting and and by setting the corresponding part of the transition matrix to , with . In our case, we used the static regression in order to prevent overfitting [79].

The parameter estimation and the model simulation are conducted in a Bayesian framework so that empirical priors can be incorporated into the model parameters [56]. Let generically denote the set of all model parameters, and let denote the full state sequence. A prior distribution on the model parameters and a distribution on the initial state values are needed to define the model. The parameters of the model in equations (C7) and (C8) are the set of variances, for which a commonly used prior distribution is the *gamma* distribution with expectation :

The hyperparameters can be interpreted as a prior sum of squares so that is a prior estimate of , and is the weight, in units of prior sample size, assigned to the prior estimate. The values of and can be sampled from using Markov Chain Monte Carlo; subsequently, the counterfactual time series are sampled from the predictive posterior distribution (see [56] for more details on the model estimation).

For this method to correctly work, the covariates themselves must not be affected by the intervention, and the relationship between covariates and treated time series, as established during the preintervention, must remain stable throughout the postintervention period. In particular, this means that the meteorological/seasonal variables used as covariates are considered to be not affected by the lockdown intervention, which seems a reasonable assumption. Given this assumption and considering that the relation between meteorological conditions and concentration remains stable after the lockdown, a significant deviation from the counterfactual concentration would indicate the causal impact of the intervention. This means that the lockdown would be considered as the event causing the drop in air pollution through activity restrictions. On the contrary, if the main drivers of the concentration were the meteorological/seasonal conditions (and not human activity), we would expect the data to follow the counterfactual prediction in the post-treatment period. In that case, no significant causal impact of the intervention would be found. Thanks to this approach, we stated the impact of the reduced human mobility and energy consumption due to the lockdown, considered as the attributable intervention on concentrations.

#### Data Availability

The data are available at https://scihub.copernicus.eu/, https://cds.climate.copernicus.eu/, and https://www.terna.it/en/electric-system/transparency-report, https://www.google.com/covid19/mobility/Accessed:<07-29-2020>, https://www.apple.com/covid19/mobility and https://data.humdata.org/dataset/movement-range-maps.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.