A Comparison Analysis of Causative Impact of PM 2.5 on Acute Exacerbation of Chronic Obstructive Pulmonary Disease (COPD) in Two Typical Cities in China

: Chronic obstructive pulmonary disease (COPD) is a major and increasingly prevalent respiratory health problem worldwide and the ﬁne particulate matter (PM 2.5 ) is now becoming a rising health threat to it. This study aims to conduct a comparison analysis of health effect on acute exacerbation of COPD (AECOPD) associated with PM 2.5 exposure in two typical cities (Beijing and Shenzhen) with different levels of PM 2.5 pollution. Both correlational relationship and causal connection between PM 2.5 exposure and AECOPD are investigated by adopting a time series analysis based on the generalized additive model (GAM) and convergent cross mapping (CCM). The results from GAM indicate that a 10 µ g/m 3 increase in PM 2.5 concentration is associated with 2.43% (95% CI, 0.50–4.39%) increase in AECOPD on Lag0-2 in Beijing, compared with 6.65% (95% CI, 2.60–10.87%) on Lag0-14 in Shenzhen. The causality detection with CCM reveals similar signiﬁcant causative impact of PM 2.5 exposure on AECOPD in both two study areas. Findings from two methods agree that PM 2.5 has non-negligible health effect on AECOPD in both two study areas, implying that air pollution can cause adverse consequences at much lower levels than common cognition. Our study highlights the adverse health effect of PM 2.5 on people with COPD after exposure to different levels of PM 2.5 and emphasizes that adverse effect in area with relative low pollution level cannot be overlooked. Governments in both high-pollution and low-pollution cities should attach importance to the adverse effects of PM 2.5 on humans and take corresponding measures to control and reduce the related losses.


Introduction
Chronic obstructive pulmonary disease (COPD) is a major and increasingly prevalent respiratory health problem in the global, with a current estimate 700 million patients and 3.2 million deaths worldwide [1,2]. So far, smoking has been recognized as the most common factor in the etiology of COPD. Other associated factors, including dusts, gases, biological and chemical exposures, have also been found in both smokers and nonsmokers [3,4]. With the deterioration of air quality in recent decades, fine particulate matter (PM 2.5 ) is now becoming a rising health threat to human respiratory system [5]. PM 2.5 is an inhalable matter which can be deposited into lower bronchus and lung parenchyma and effect the respiratory system by inducing inflammatory reaction, which is known to cause acute exacerbations of respiratory disease [6]. Recently, PM 2.5 has been proven to be a primary detrimental factor for COPD by plentiful epidemiological studies [7][8][9][10]. In addition, there are also substantial evidences indicating that short-term exposure to PM 2.5 can lead to acute exacerbation and even an increased risk of death in COPD patients [11][12][13][14]. A panel study in Beijing [15] also finds out that COPD patients are more susceptible to respiratory inflammation following PM 2.5 exposure than individuals without COPD.
It is reported that the morbidity and mortality rate of COPD is increasing in the Asia-Pacific regions [16,17]. In China, the prevalence of COPD in the population over 40 years old is estimated to be 13.7% [18], making it a health issue that should not be overlooked. In addition, COPD usually brings a substantial economic burden to healthcare services [19]. According to the National Heart Lung Blood Institute, the direct cost of COPD was estimated to be $29.5 billion in the U.S. in 2009 [20]. In Singapore, the average direct cost of COPD was $9.9 million per year from 2005 to 2009 [17]. Therefore, controlling the burden of COPD can be beneficial to effectively alleviate the pressure on healthcare services and social economy. To reduce the disease burden of COPD, it is important to explore the adverse health effects on COPD patients associated with exposure to PM 2.5 pollution. Several time series studies have been conducted in this purpose. A case crossover study in Ontario, U.S. [21] observed positive associations of COPD with PM 2.5 , indicating a 7.00% (95% confidence interval, CI, 6.00-8.00%) increase in COPD exacerbation cases associated to 3.4 µg/m 3 increase in PM 2.5 concentration. Similar results were also detected in China. Ko et al. [22] observed that a 10 µg /m 3 increase in the 5-day cumulative average concentration of PM 2.5 was associated with 3.10% (95% CI, 2.60-3.60%) increases in COPD hospital visits, in Hongkong. Bao et al. [23] found that a 10 µg/m 3 increase in PM 2.5 concentration at a lag of 0-7 days was associated with 1.190% (95% CI, 0.176-2.215%) increase in hospital outpatients for COPD in Lanzhou. A newly conducted case-crossover study in Taiwan [24] showed that each interquartile range increase in PM 2.5 was associated with 6.6% (95% CI, 0.5-13.0%) increase in risk of COPD exacerbation.
However, there are still some insufficiencies in current studies. On one hand, related studies in Chinese cities are limited and lack comparison among cities with different PM 2.5 levels. On the other hand, most of the studies mainly focus on the correlational relationship between the disease and air pollutions revealed by relative risk (RR) or odds ratio, the causal connection between the two is unknown. In this study, we intend to explore the adverse health effect of PM 2.5 on acute exacerbation of COPD (AECOPD) in Beijing and Shenzhen, two typical cities with large differences in PM 2.5 concentration. In addition, to reveal the correlational relationship between AECOPD and PM 2.5 by adopting a regression model-based lag effect analysis, we also detect the causal connection between the two by applying a model-free method.

Study Areas
To perform a comparison analysis on health effect associated with different levels of PM 2.5 , two typical study areas, Beijing urban area and Shenzhen City are selected. Beijing, as the capital city of China, is the national political and economic center. With the development of Jing-Jin-Ji metropolitan area, which is dominated by heavy industry, air quality in Beijing has suffered a severe deterioration in recent decades [25]. Haze weather occurs frequently and the adverse health effect of air pollution, especially PM 2.5 , has drawn great public attention. Therefore, Beijing urban area (including Dongcheng, Xicheng, Chaoyang, Haidian, Fengtai, Shijingshan), is chosen as the typical area with high PM 2.5 level in this study for comparison. Shenzhen, a young city in China, was once praised as an Environmental Protection Model City in 1997 for its favorable air environment. Playing an important role in the Pearl River Delta region, which is one of the most developed regions with the highest aggregation of industry in China, air quality in Shenzhen has also deteriorated and affected the living condition of local citizens in some way [26,27].
Nevertheless, compared to Beijing, PM 2.5 level in Shenzhen is still considered quite low [28]. In this study, Shenzhen City is chosen as the typical area with low PM 2.5 for comparison.

Data Collection
Daily AECOPD records were obtained at ten comprehensive hospitals in Beijing urban area (from Xu's study [29]) and 66 major hospitals in Shenzhen City (Figure 1), from 1 January 2013 to 31 December 2013. Monitoring data of PM 2.5 was gathered from Beijing Environmental Protection Bureau and Shenzhen Environmental Monitoring Center, including data from 17 monitoring sites located in Beijing urban area and 19 sites in Shenzhen City (Figure 1), following the Chinese National Ambient Air Quality Standards (GB3095-2012) [30]. Daily-averaged concentration is used to represent the degree of the exposure to PM 2.5 . There are some missing values in the 17 monitoring sites in Beijing which may cause some gaps in the time series analysis. To fill them up, we applied the linear interpolation method by using a three-day average. In addition, to control for the confounding factors, we also collected the weather conditions during the same period, involving daily mean temperature ( • C) and relative humidity (%), from the official website of the Chinese Meteorological Bureau.
this study for comparison. Shenzhen, a young city in China, was once praised as an Environmental Protection Model City in 1997 for its favorable air environment. Playing an important role in the Pearl River Delta region, which is one of the most developed regions with the highest aggregation of industry in China, air quality in Shenzhen has also deteriorated and affected the living condition of local citizens in some way [26,27]. Nevertheless, compared to Beijing, PM2.5 level in Shenzhen is still considered quite low [28]. In this study, Shenzhen City is chosen as the typical area with low PM2.5 for comparison.

Data Collection
Daily AECOPD records were obtained at ten comprehensive hospitals in Beijing urban area (from Xu's study [29]) and 66 major hospitals in Shenzhen City (Figure 1), from 1 January 2013 to 31 December 2013. Monitoring data of PM2.5 was gathered from Beijing Environmental Protection Bureau and Shenzhen Environmental Monitoring Center, including data from 17 monitoring sites located in Beijing urban area and 19 sites in Shenzhen City (Figure 1), following the Chinese National Ambient Air Quality Standards (GB3095-2012) [30]. Daily-averaged concentration is used to represent the degree of the exposure to PM2.5. There are some missing values in the 17 monitoring sites in Beijing which may cause some gaps in the time series analysis. To fill them up, we applied the linear interpolation method by using a three-day average. In addition, to control for the confounding factors, we also collected the weather conditions during the same period, involving daily mean temperature (°C) and relative humidity (%), from the official website of the Chinese Meteorological Bureau.

Lag Effect Analysis
The core model used in lag effect analysis is a generalized additive model (GAM) which is suitable for exploring and fitting non-linear relationship between the response and the predictor variables. In this study, we use it to investigate the exposure-response relationship between the exposure to PM 2.5 and the occurrence of AECOPD. However, the adverse health effect of the exposure is not limited to the period when the exposure occurs. There could be a lag before a response in human body. Therefore, we combined the GAM with a distributed lag nonlinear model (DLNM) [31] which can provide a detailed time-course representation of the lag procedure, to examine the health effect associated with PM 2.5 . The main advantage of the DLNM is that it allows the model to estimate the overall cumulative lag effect as the sum of the single-day lag effects over the whole selected lag period [32]. In addition, the level of PM 2.5 , we also added several confounding variables involving the calendar time, week and holiday effect, as well as meteorological variables for improving the model performance [33,34]. The calendar time, temperature and relative humidity are expressed as smoothing spline functions, while the week and holiday effects are added as two dummy variables in the model. Degrees of freedom (df ) of smoothing functions are determined by the Akaike's information criterion. We applied a generalized cross validation to guide the determination of df until the sum of absolute difference got a minimum. In this study, df was set to 7 for calendar time and 3 for temperature and relative humidity in both two study areas to account for the potential nonlinear effects. The final lag effect model is of the following form: where t is the day for observation; X t−L represents the concentration of PM 2.5 with a lag period of L; f (X t−L ) stands for the lag model involving a single-day lag mode and a cumulative lag mode; E(Y t ) is the expected number of cases on day t; α is the intercept term; β represents the log-RR of cases associated with a unit increase of PM 2.5. ; S(time, d f ) and S(Z t , d f ) are the smoothing splines of calendar time and metrological factors (temperature and relative humidity); DOW and Holiday stand for the week and holiday effects with coefficient β w and β h , respectively.
To perform the analysis, DLNM and Mixed GAM Computation Vehicle (MGCV) packages in R (3.4.0, University of Auckland, Auckland, New Zealand) were used. The outcomes are calculated as RR or percent change in daily counts and at the 95% CI, in association with a unit (10 µg/m 3 ) increase of PM 2.5 concentration.

Causative Impact Detection
The causative impact detection is conducted by applying the convergent cross mapping (CCM) method proposed by Sugihara et al. [35]. Previous studies prove that it shows advantage in distinguishing causality from standard correlations based on a time series analysis (at least 25 pairs of observations) [36,37]. In nonlinear dynamic systems, variables are often inseparable and contain multi-dimensional information, making unstable correlation and weak coupling relation common in the correlation relationships of the system. For the most frequently used causality analysis method, the Granger causality analysis [38], it tends to behave poorly in detecting causality in weak-to-moderate coupling relations. Instead, CCM is insensitive to manual setting of parameters and can extract reliable causality between diverse variables, making it more suitable for revealing potential causality in nonlinear dynamic systems [39,40]. In CCM, the complex and nonlinear systems is analyzed through state-space reconstruction and the reconstructed shadow manifold of the two variables corresponds to each other in the time dimension. Giving two time series of variables X{x 1 ,x 2 , . . . ,x L } and Y{y 1 ,y 2 , . . . ,y L } (L is the length of time period), CCM algorithm can be conducted in four steps, as follows: Step One Reconstruct the shadow manifold, M x and M y , using the lagged-coordinate vectors X and Y, which is: Atmosphere 2021, 12, 970 5 of 12 where t ranges from 1 + (E − 1) τ to L; E is the dimension in which the reconstructed attractor is embedded; τ is the time lag. Assuming there is N pairs of observations, E ≤ N − 1 [35].
Step Two Determine the weight. Since M y is diffeomorphic to M x , there will be synchronous lagged-coordinate vector x(t) and its E + 1 nearest neighbors on M x , which can be used to build the weight w i , defined as: where For each x(t), the nearest neighbor search gets a set of distances sorted from the closest to the outermost by an associated set of time {t 1 ,t 2 , . . . ,t E+1 }. The distance is measured as the Euclidean distance between the two vectors.
Step Three Use the w i to create a cross-mapped estimation of y t by calculating with a weighted mean the nearest neighbors in M y . The cross-mapped estimation is express as:ŷ Step Four Calculate the CCM correlation coefficient. The degree of consistency between the cross-mapped estimation and the true value determines the predictive ability of Y on X [37], which can be quantified by Pearson correlation coefficient between original and estimated time-series. It is also referred to as the CCM correlation coefficient (ρ), defined as: Meanwhile, a t-statistic for correlation coefficient at a level of significance is calculated as: In CCM method, determining the optimal E and τ is crucial to the analysis. Supposing that E max is the optimal dimension, the dimensionality generically ranges between (E max − 1)/2 and E max , based on Whitney's theorem [41]. In this study, dimension (E) from the two times series is defined as 2 using the false nearest neighbor method and the value of τ is set to 2 based on the average mutual information criterion for both study areas.

Descriptive Statistics
During the study period, the average daily PM 2.5 concentration was 102 µg/m 3 , ranging from 6.7 µg/m 3 to 508.5 µg/m 3 in Beijing and 37 µg/m 3 , ranging from 7.9 µg/m 3 to 129.  Table 1. The average daily count of AECOPD was 2 in Beijing (ranged from 0 to 7) and 10 in Shenzhen (ranged from 1 to 23). The time series graph in Figure 2 shows the temporal distribution of AECOPD and PM 2.5 concentrations in 2013. The curve of PM 2.5 in Shenzhen showed an apparent decrease of concentration during the summer period (June, July, August), which is correlated to the monsoon from the sea in summer [42]. The curve of PM 2.5 in Beijing showed no noticeable trend. Judging from the time series graph, no particular association between the AECOPD and PM 2.5 can be observed in both study areas through the unaided viewing. Shenzhen when the daily PM2.5 concentration met the national standard, according to the Chinese Ambient Air Quality Standards (Grade II, 75 μg/m 3 for 24 h average PM2.5 concentration). However, only 30 days in Beijing and 120 days in Shenzhen met the WHO Air Quality Standards, which is 25 μg/m 3 for 24 h average PM2.5 concentration. Based on the statistics above, PM2.5 pollution in Beijing was much more severe than that in Shenzhen in 2013. At the same time, a total of 712 and 3634 AECOPD cases were recorded from the selected hospitals in Beijing and Shenzhen in 2013, respectively. The statistical characteristics of AECOPD, PM2.5 concentration and meteorological factors in Beijing and Shenzhen in 2013 are summarized in Table 1. The average daily count of AECOPD was 2 in Beijing (ranged from 0 to 7) and 10 in Shenzhen (ranged from 1 to 23). The time series graph in Figure 2 shows the temporal distribution of AECOPD and PM2.5 concentrations in 2013. The curve of PM2.5 in Shenzhen showed an apparent decrease of concentration during the summer period (June, July, August), which is correlated to the monsoon from the sea in summer [42]. The curve of PM2.5 in Beijing showed no noticeable trend. Judging from the time series graph, no particular association between the AECOPD and PM2.5 can be observed in both study areas through the unaided viewing.

The Lag Effect of PM 2.5 on AECOPD
The lagged health effect on AECOPD associated with short-term exposure to PM 2.5 was analyzed using the GAM model introduced in Section 2.3.1. The lag period was set to range from 0 to 14 days (denoted as Lag0, Lag1 . . . Lag14) in both two study areas. In the model, effects of PM 2.5 on AECOPD was considered as a linear correlation relationship. The RR for AECOPD was calculated as variation in counts of daily cases associated with 10 µg/m 3 increase of PM 2.5 concentration. The exposure-response graphs for single-day lag effect and cumulative lag effect are shown in Figure 3, in which the 95% CIs are represented as black bars and gray areas. The cumulative lag effect was computed by applying a polynomial curve fitting on single day lag effect at a degree of 3 in both study areas.

The Lag Effect of PM2.5 on AECOPD
The lagged health effect on AECOPD associated with short-term exposure to PM2.5 was analyzed using the GAM model introduced in Section 2.3.1. The lag period was set to range from 0 to 14 days (denoted as Lag0, Lag1 … Lag14) in both two study areas. In the model, effects of PM2.5 on AECOPD was considered as a linear correlation relationship. The RR for AECOPD was calculated as variation in counts of daily cases associated with 10 μg/m 3 increase of PM2.5 concentration. The exposure-response graphs for single-day lag effect and cumulative lag effect are shown in Figure 3, in which the 95% CIs are represented as black bars and gray areas. The cumulative lag effect was computed by applying a polynomial curve fitting on single day lag effect at a degree of 3 in both study areas.   Tables 2 and 3 show the specific RR values for AECOPD associated with a 10 μg/m 3 increase in PM2.5 in terms of single-day lag effect and cumulative lag effect, together with the 95% CIs in the two study areas. For Beijing urban area, significant singleday lag effects are only observed at Lag0 (the same day) and Lag1, with a maximum RR of 1.0155 (95% CI, 1.0058-1.0253) at Lag0; and significant cumulative lag effects are observed from Lag0 to Lag3, with a maximum cumulative RR of 1.0243 (95%, 1.0050-1.0439) at Lag2. For Shenzhen city, significant single-day lag effect appears from Lag2 to Lag8, with a maximum RR of 1.0065 (95% CI, 1.0010-1.0120) at Lag5; and significant cumulative  Tables 2 and 3 show the specific RR values for AECOPD associated with a 10 µg/m 3 increase in PM 2.5 in terms of single-day lag effect and cumulative lag effect, together with the 95% CIs in the two study areas. For Beijing urban area, significant singleday lag effects are only observed at Lag0 (the same day) and Lag1, with a maximum RR of 1.0155 (95% CI, 1.0058-1.0253) at Lag0; and significant cumulative lag effects are observed from Lag0 to Lag3, with a maximum cumulative RR of 1.0243 (95%, 1.0050-1.0439) at Lag2. For Shenzhen city, significant single-day lag effect appears from Lag2 to Lag8, with a maximum RR of 1.0065 (95% CI, 1.0010-1.0120) at Lag5; and significant cumulative lag effect appears from Lag4 and remains significant until Lag14, with a maximum cumulative RR of 1.0665 (95% CI, 1.0260-1.1087) at Lag14. Though significant lag effects on AECOPD associated with increase in PM 2.5 concentration are detected in both two study areas, the lag patterns are quite different. The response time of the lag effect in Beijing is shorter than that in Shenzhen, which occurs on the same day of exposure and lasts for a small period (3 days in cumulative mode). In contrast, the lag effect in Shenzhen does not emerge until two days after the exposure and the cumulative lag effect could last for at least 10 days. Judging by the value of single-day RR, the same amount of increase in PM 2.5 concentration has slightly greater effect on AE-COPD in Beijing urban area than in Shenzhen. However, for the cumulative lag effect, the same amount of increase in PM 2.5 concentration is associated with much more responses on AECOPD in Shenzhen than in Beijing urban area.

Results in
To sum up, PM 2.5 exposure was proved to have non-negligible lag effect on AECOPD in both two study areas, regardless of the level of PM 2.5 concentration. The regression model-based method applied in this section reveals the correlational relationship between the two variables. However, correlation does not mean causality [39] and the causal connection between AECOPD and PM 2.5 exposure still remained unconfirmed. In the following section, we will discuss the causative impact of PM 2.5 on AECOPD based on the results from a model-free analysis by using CCM method.

The Causative Impact of PM 2.5 on AECOPD
To confirm the causality using CCM method, two conditions are required: (1) the convergence result of the CCM correlation coefficient (ρ) is significantly larger than 0; (2) the pattern of CCM curve shows an obvious increasing and convergence trend with the growth of the time series length (L). The convergent cross maps of Beijing urban area and Shenzhen City are shown in Figure 4, in which red lines represent causative impact of AECOPD count on PM2.5 concentration while the blue lines represent the opposite.

The Causative Impact of PM2.5 on AECOPD
To confirm the causality using CCM method, two conditions are required: (1) the convergence result of the CCM correlation coefficient ( ) is significantly larger than 0; (2) the pattern of CCM curve shows an obvious increasing and convergence trend with the growth of the time series length (L). The convergent cross maps of Beijing urban area and Shenzhen City are shown in Figure 4, in which red lines represent causative impact of AECOPD count on PM2.5 concentration while the blue lines represent the opposite.
As shown in Figure 4, the blue lines in the two graphs meet the required conditions of CCM, indicating that PM2.5 concentration has causative impacts on the counts of AECOPD in both Beijing urban area and Shenzhen City. Regarding the red lines in two graphs, the required conditions are not fulfilled, implying that count of AECOPD has no causal impact on PM2.5 concentration. Numerical results of causality ( value) in Table 4 illustrate the strength of the causative impact and similar causality strength is observed for both study areas.   In general, significant causative impact of PM2.5 exposure on AECOPD has been confirmed in both two study areas, regardless of the level of PM2.5 concentration. Different from the common cognition that exposure to higher level of PM2.5 may have greater effect on human health, this study finds that the effect of PM2.5 on human body has little to do with the concentration level. In most cases, people are likely to underestimate the adverse health effect of PM2.5 when the exposure level is low. However, this study provides a persuasive evidence to show that health effect in area with relative low pollution level also need to be taken seriously, especially for people with chronic respiratory disease. As shown in Figure 4, the blue lines in the two graphs meet the required conditions of CCM, indicating that PM 2.5 concentration has causative impacts on the counts of AECOPD in both Beijing urban area and Shenzhen City. Regarding the red lines in two graphs, the required conditions are not fulfilled, implying that count of AECOPD has no causal impact on PM 2.5 concentration. Numerical results of causality (ρ value) in Table 4 illustrate the strength of the causative impact and similar causality strength is observed for both study areas. Table 4. Numerical result of causality (ρ value) between PM 2.5 and AECOPD in Beijing urban areas and Shenzhen City. (Value with * is significant at the 0.05 level).

Combination
Causative Impact (ρ Value) In general, significant causative impact of PM 2.5 exposure on AECOPD has been confirmed in both two study areas, regardless of the level of PM 2.5 concentration. Different from the common cognition that exposure to higher level of PM 2.5 may have greater effect on human health, this study finds that the effect of PM 2.5 on human body has little to do with the concentration level. In most cases, people are likely to underestimate the adverse health effect of PM 2.5 when the exposure level is low. However, this study provides a persuasive evidence to show that health effect in area with relative low pollution level also need to be taken seriously, especially for people with chronic respiratory disease.

Conclusions
In this study, a comparison analysis of health effect associated with PM 2.5 is conducted in two typical cities with different levels of PM 2.5 concentration. Both correlational relationship and causal connection between PM 2.5 exposure and AECOPD are investigated by adopting a time series analysis based on GAM and CCM. Results from the two methods agree that PM 2.5 exposure has non-negligible health effect on AECOPD in both two study areas, implying that air pollution cause adverse consequences at much lower levels than common cognition. In conclusion, our study highlights the adverse health effects of PM 2.5 pollution on people with chronic respiratory disease in cities with both high and low levels of PM 2.5 in China and emphasizes adverse effect in area with relative low pollution level cannot be overlooked. Governments in both high-pollution and low-pollution cities should attach importance to the adverse effects of PM 2.5 on human body and take corresponding measures to control and reduce the related losses.