Abstract
Objective
Hemorrhagic fever with renal syndrome (HFRS), a natural–focal infectious disease caused by hantaviruses, resulted in 37 deaths between 2011 and 2015 in Hubei Province, China. HFRS outbreaks are seasonally distributed, exhibiting heterogeneity in space and time. We aimed to identify the spatial and temporal characteristics of HFRS epidemics and their probable influencing factors.
Methods
We used the space–time cube (STC) method to investigate HFRS epidemics in different space–time locations. STC can be used to visualize the trajectories of moving objects (or changing tendencies) in space and time in three dimensions. We applied space–time statistical methods, including space–time hot spot and space–time local outlier analyses, based on a calculated STC model of HFRS cases, to identify spatial and temporal hotspots and outlier distributions. We used the space–time gravity center method to reveal associations between possible factors and HFRS epidemics.
Results
In this research, HFRS cases for each space–time location were defined by the STC model, which can present the dynamic characteristics of HFRS epidemics. The STC model delivered accurate and detailed results for the spatiotemporal patterns of HFRS epidemics.
Conclusion
The methods in this paper can potentially be applied for infectious diseases with similar spatial and temporal patterns.
Keywords
Introduction
Hemorrhagic fever with renal syndrome (HFRS), characterized by fever, hemorrhage, kidney damage, and hypotension, is an important infectious disease that is mainly caused by Hantaan orthohantavirus (HTNV) and Seoul orthohantavirus (SEOV).1–4 HFRS is distributed in several countries, with China the most seriously affected. HFRS in China has accounted for more than 90% of cases worldwide, 5–9 and Hubei Province has the highest occurrence. Since the first case of HFRS emerged in 1957 in Hubei, HFRS epidemics have expanded, reaching a peak in 1983 with 23,943 cases. During the period 1980–2009, a total 104,467 cases of HFRS have been reported in Hubei Province. The increase in HFRS cases has a considerable impact on human health and social stability.10–12
Spatial and temporal statistical methods have been used worldwide to determine the spatiotemporal distributions and clustering characteristics of HFRS,12,13 such as in Buenos Aires Province of Argentina, 14 Germany, 15 the city of Brussels, and across Belgium. 16 In China, the Moran’s I index has been adopted in spatial global autocorrelation analysis to identify the overall spatiotemporal pattern of HFRS outbreaks in Hubei Province, China during 2005–2014. 17 The Kulldorff spatial scan statistic has been used to identify clustering of HFRS using data during 1980–2009. 18 These methods enable interactive examination of global and local clustering patterns in a specific temporal profile; however, connections and influences from different profiles are not considered. Liu et al. 19 applied the ARIMA (Autoregressive Integrated Moving Average) model to evaluate and forecast the HFRS incidence in China using time series historical HFRS incidence data. The ARIMA model provides an important tool for estimation using continuous time series data; however, the correlations for different spatial locations are ignored in this analysis. To address this limitation, Martin and Oeppen (1975) proposed the STARMA (Space-Time Autoregressive Moving Average) model and extensions for spatial and temporal data. 20 Nonetheless, the spatial influences of data in distinct locations are only considered as individual factors in the STARMA model. The accuracy of analysis for spatial epidemics like of HFRS may be further improved by comprehensive consideration of both spatial and temporal characteristics, as well as interactive influences.
We created a three-dimensional map of HFRS cases to visually and simultaneously comprehend both the temporal duration and spatial extent of HFRS cases in each space–time location. This was accomplished using a space–time cube (STC) method. Spatial and temporal clustering methods based on the STC model, which has the advantage of tracking the spatial and temporal trajectories of moving objects or tendencies, have been used to identify the dynamic tendencies and patterns of HFRS cases. The STC was introduced by Hagerstrand. Rucker, and Szego 21 who proposed the method to create a three-dimensional cube that can contain spatial and temporal data for specified data requiring further space–time pattern mining. The STC method has evolved since its introduction; it has also been used to examine time series trends across a studied area. 22 Recently, the STC model has been applied to reveal the spatiotemporal patterns in a data set for users in different areas and events. 23 This method can provide important support for studies in epidemiology, health geography, and environmental regulation. 24 For instance, the STC method combined with space–time scan statistics was applied to analyze and display the spatiotemporal pattern of incidence of hand-foot-and-mouth disease (HFMD) in Guangdong Province from May 2008 to March 2009. 25
In the present study, we first applied the STC and its clustering methods to detect spatial and temporal dynamic patterns of HFRS cases in Hubei Province. The analysis of HFRS in this work consisted of three main steps. In the first step, we adopted STC-based approaches, including emerging hot spot analysis and space–time local outlier analysis, to investigate the spatial and temporal clustering characteristics and locations of HFRS cases. The second step focused on clarifying the key influencing factors of the disease in various regions. The third step used a spatial gravity center model and correlation analysis to analyze the moving trajectory of HFRS. Last, possible influencing factors of HFRS epidemics were additionally explored.
Study data and methodology
Study data
The study data consisted of geographic, health (i.e., number of HFRS cases), climate, and human population density data. These data span a 5-year period (2011–2015) in Hubei Province, China. Geographic information data were obtained from the China National Administration of Surveying, Mapping and Geoinformation. The Hubei Province Center for Disease Control and Prevention of the Chinese Center for Disease Control and Prevention was the source for HFRS data (number of cases per month in the county). Climate data were obtained from the National Centers for Environmental Prediction and Hubei Meteorological Bureau. The Hubei Statistical Yearbook was the source of annual population data. Esri ArcMap 10.5 was used to organize and process the geographic data.
Methodology
Figure 1 depicts the relationships between the statistics used in this study.

Relationships between the statistics used in this study Figure represents the overall analysis process.
Construction of space–time cube (STC)
The STC method has been used to create a three-dimensional cube containing spatial and temporal data for a specified data set requiring further space–time pattern mining. In the present study, we used the STC model to aggregate HFRS cases into a series of spatiotemporal bins (refer to Figure 2). Each bin is characterized in space, using a fixed (x,y) coordinate, and time (t). Bins that cover the same area share a common location ID; bins that share the same duration have a common time-step ID. We used different bins to represent the number of HFRS cases in different space–time locations. Data in a row of the cube collection represent the case values in different locations with a specific date; data in a column represent the time series for a specific space location. Continuous statistical methods, including emerging hot spot analysis and local outlier analysis, were applied to examine the spatial and temporal dynamic patterns of HFRS epidemics.

Space–time cube graphic sample Figure gives an explanation for each value in the space–time cube.

netCDF Classic Data Model Figure depicts the data structure for the netCDF data format.
The steps of the STC method in the analysis are as follows:
Aggregate the collections of points into space–time bins and integrate the bins into a netCDF data model (refer to Figure 3). The netCDF data set includes attributes, variables, and dimensions. The bins are uniquely identified using a name and ID number.
26
Calculate each point and aggregate the attributes of HFRS cases. Perform the Mann–Kendall (M-K) trend test for every county with data. Each of these is treated as an independent bin time series test.
In this work, the M-K trend test is applied independently to each location that has data. The data for each location are stored in a bin with a value (i.e., count) for one period of time; a series of values for the bin is recorded over time. The computation of the M-K trend test uses straightforward, pairwise comparisons of bin values over time. If the bin value for a time period is smaller than that for the next time period, then the sum is incremented. If the result is larger, then the sum is decremented. If the values are the same (a tie), the sum remains unchanged. After the pairwise comparisons are complete, the expected sum is 0, which indicates that no positive or negative trend exists in the data. Actual trends in the data are computed using the variance, number of ties and time periods, (observed) sum, and the expected sum. The z-score and p-value are computed for each time series. The sign of the z-score indicates the positive or negative direction of the trend; 27 the p-value indicates its statistical significance.
Based on the calculated STC results, we conducted space–time hot spot analysis and space–time local outlier analysis to further explore the spatial and temporal tendencies of HFRS in Hubei Province. In this study, emerging hot spot analysis was used to identify the trend in the spread of HFRS cases; local outlier analysis was used to identify significant clusters and outliers.
Spatial and temporal weight matrix
For many statistics, including the Getis–Ord Gi* and Anselin local outlier analysis, the spatial relationships between each location are specified through a spatial weight matrix. However, when conducting spatial and temporal-based statistical analysis, the spatial weight matrix is insufficient to express the spatial and temporal relationships for each location. In this case, the spatial and temporal weight matrix is extended. There are commonly three types of spatial and temporal weight matrices. First, based on the rules of the spatial weight matrix, a specified spatial distance d and a specified time interval t are set up. Only when there is a spatial distance smaller than d and temporal interval smaller than t, elements in the matrix are set to 1 and all others are 0 (Smith and Wu, 2009). Second, the spatial weight matrix is extended, a spatial and temporal distance proposed, and values are calculated for the matrix based on the spatial and temporal distance in a specified threshold. Third, an external temporal matrix is built and the spatial and temporal weight matrix with the spatial matrix set up through a Kronecker product. In this study, the spatial and temporal distance is determined using the second mode mentioned above. The spatial and temporal weight matrix is calculated using Equation 1.
In Equation 1,
In Equation 3,
Spatial and temporal-based Gi* statistic
In this work, we proposed the spatial and temporal based Getis–Ord Gi* statistic (ST-G*). Spatial and temporal weight is calculated to obtain the Getis–Ord Gi* statistic. Compared with traditional spatial–temporal distribution analysis, analysis using the ST-G* has the advantage of analyzing continuous spatial and temporal trends. With the created STC, the ST-G* can calculate the Getis–Ord Gi* statistic with the provided neighborhood distance and neighborhood time-step parameters. In addition to existing hot spots, the results of calculation reveal new hot spots that are changing (e.g., intensifying or diminishing) or fluctuating between hot and cold. We used the ST-G* to identify trends in the clustering of HFRS cases for each county location, based on the STC model. To carry out ST-G* analysis, the ST-G* statistics of each bin were computed based on the spatial and temporal weight matrix. The ST-G* statistic was used to explore the local spatial autocorrelation of the collected data sets and to detect the spatial differences caused by spatial autocorrelation. From these results, hot and cold spot areas of the spatial object attribute distribution can be judged. 12
The ST-G* statistic introduces the concept of spatiotemporal proximity; the formulas are as follows7,28
The result of regression is directly impacted by the computed value of the bandwidth. We used the Akaike information criterion (AIC) to identify the bandwidth for each space–time bin. 29 The Gi* calculated with Equation 4 is the z-score. For statistically significant positive z-scores, the larger the z-score, the more intense the clustering of high values (i.e., a hot spot).30,31
Spatial and temporal-based Anselin local Moran's I
We used the spatial and temporal-based Anselin local Moran's I (ST-ALM) to identify statistically significant clusters and outliers in the context of both space and time for a study area. This analysis identified the locations exhibiting distinct statistical differences in comparison with their neighbors.
Based on the STC model, the ST-ALM statistic for each bin in the cube was computed. For space–time bin
Space gravity center model
To explore the correlation between the influencing factors of population migration and HFRS epidemics, we also investigated the spatial gravity center of HFRS cases and spatial gravity center of the population distribution in this study. The space gravity center coordinate is calculated as follows35,36
Results and analysis
Features of HFRS outbreaks in Hubei Province
Figure 4 demonstrates the time series of HFRS incidence in Hubei Province from 2011 to 2015. The figure shows that the HFRS incidence fluctuated and had no decreasing trend during this period. Obvious peaks and troughs can be detected in the figure, which means that HFRS epidemics had a significant seasonal distribution pattern. A nonlinear distribution pattern appears in the curve, from which it can be inferred that HFRS incidence from 2011 to 2015 had no statistically significant temporal distribution.

HFRS incidence in Hubei Province from 2011 to 2015 Each dot in the figure represents the monthly provincial incidence of HFRS epidemics.
The spatial distribution of annual HFRS cases is presented in Table 1. In four of the five years, the HFRS epidemic presented a clustered pattern of distribution in space, with a z-score >1.65; these results were statistically significant, with p < 0.05. The results indicate that the occurrence of HFRS in Hubei Province is globally autocorrelated and temporally continuous. These results support the execution of STC-based analysis.
Spatial autocorrelation results of annual HFRS cases (Hubei Province, 2011–2015).
HFRS, hemorrhagic fever with renal syndrome.
STC-based analysis
Based on the above characteristics of cases, we developed the STC model of HFRS cases with a time step of 1 month and space distance of 10,000 meters. This model contains 60 time-step intervals, and the number of space–time bins is 4,000. The entire length of time covered by the HFRS cases is 60 months (January 2011 to December 2015). The M-K statistic test determines the final bin counts. Any bins that cannot be filled based on the estimation criteria result in the whole location being excluded from the analysis. The calculating attributes used is HFRS cases; the statistic fill method is space–time neighbors, which means empty bins are filled using a interpolated univariate spline algorithm. The M-K statistic test resulted in 3,300 bins (82.5%) included in the analysis.
After constructing the STC model, the number of HFRS cases at each position in every time-step interval was set as the time series. We applied M-K statistical trend analysis; this technique is based on a spatial time series analysis. It analyzes the number of points for all locations in each time-step interval as a time series of count values, based on the M–K statistic. The z-score was –5.709 and the corresponding p-value was 0. This result suggests that Hubei Province shows a statistically significant decreasing trend in the overall number of HFRS cases.
Results of ST-G*
Figure 5 shows the emerging hot spot distribution and the statistical results. According to Figure 5, the following conclusions can be drawn.

Visualization of the results of hot spot analysis: HFRS cases in Hubei Province (2011–2015) Space–time cubes represented by the analyzed spatial and temporal patterns.
Zhongxiang is a diminishing hot spot county. This category describes a statistically significant hot spot for 90% or more of the time-step intervals, including the final time step, with a significantly decreasing change over time. For this type of county, the number of HFRS cases is temporally increasing, with the incidence risks ranking the highest of all hot spot types. Thus, prevention of HFRS in this county should receive greater attention, to slow the spread of HFRS epidemics.
Yicheng and Zaoyang are persistent hot spot counties. This category describes a statistically significant hot spot for 90% or more of the time-step intervals, with no significant change over time. For these two counties, the number of HFRS cases temporally remains at a high level without an increasing tendency. Thus, prevention and treatment of HFRS in these counties could contribute to reducing the overall HFRS incidence.
Zuogui, Changyang, and Fangxian are intensifying cold spot counties, which means a statistically significant cold spot for 90% or more of the time-step intervals, including the final time step. In addition, the clustering intensity of the low counts in each time step is increasing significantly. For these counties, the number of HFRS cases is decreasing gradually; hence, HFRS epidemics have been effectively controlled in these counties.
Shiyan, Guangshui, Xiaogan, and Badong are sporadic cold spot counties, which are intermittent cold spots. These are statistically significant cold spots for less than 90% of the time-step intervals and none of the time-step intervals have been statistically significant hot spots. For these counties, HFRS epidemic risks still exist and the counties cannot be regarded as noninfected areas.
The M-K trend test was performed on every location with HFRS case data as an independent bin time series test. Locations with bins indicated that HFRS cases in the location showed significant temporal trends. However, 45 of the 55 calculated STC models showed no detectable pattern. For these counties, other data analysis methods should be used to determine the spreading tendencies of HFRS epidemics. The ST-G* analysis results indicated that the epidemic hot spots are clustered in the central part of Hubei Province, which means there may be favorable environments and risk factors for pathogens to grow and reproduce in these counties.
Results of ST-ALM
Local outlier analysis demonstrated that 34 of 55 locations have outliers. Based on the results in Figure 6 and Table 2, the following conclusions can be drawn.

Visualization of local outlier analysis results: HFRS cases in Hubei Province (2011–2015) Space–time cubes represented by the analyzed spatial and temporal outliers.
Local outlier analysis summary results.
Five counties (Yunxian, Zhuxi, Anlu, Yingcheng, and Jingshan) showed no significant clusters or outliers, indicating spatially and temporally random distributions of HFRS epidemics. No counties showed a high HFRS incidence over time; 156 bins demonstrated only high-high clusters, and no counties had only statistically high-high clusters. Four counties (Gucheng, Laohe, Tianmen, and Xiantao) showed only low-high outliers, indicating that these are low-incidence counties surrounded by high-incidence counties. However, prevention of HFRS is still needed, with qualified monitoring.
Sixteen counties (Zhushan, Shiyan, Danjiang, Fangxian, Baokang, Shennong, Xingshan, Zuogui, Changyang, Yichang, Changyang, Zhijiang, Yidu, Wufeng, Guangshui, and Hong’an) demonstrated only low-low clusters. These counties can therefore be regarded as having a low risk of HFRS epidemics. Six counties (Yunxi, Suizhou, Maocheng, Yicheng, Lichuan, and Enshi) showed only high-low outliers, indicating that these are surrounded by counties with a low HFRS incidence. Therefore, more attention is needed for counties near these high-low counties, for early prevention of HFRS.
A total 24 counties were categorized as having multiple types of statistically significant clusters and outliers over time. Therefore, further analysis should be conducted for these counties. From this analysis, it can be inferred that no statistically high-high clusters exist, which indicates that HFRS epidemics in Hubei Province are effectively controlled.
Discussion
In this study, we investigated the spatiotemporal distribution and clustering patterns of HFRS cases in Hubei Province. We carried out an analysis of the correlation between HFRS cases and possible influencing factors including landscape characteristics, population, average humidity, rainfall, and average temperature. This analysis was conducted according to two methods: linear regression analysis using the ordinary least square (OLS) method, and nonlinear regression analysis of the geography using geographically weighted regression (GWR).37,38
As shown in Table 3, the OLS analysis demonstrated that farmland area, water area, population, average humidity, and rainfall are positively correlated with HFRS cases whereas forest area and average temperature had no significant correlation. Compared with the results demonstrated using the OLS method, the results from GWR analysis indicated that forest area and average temperature are also related to the emergence of HFRS whereas population factors had no significant effect on the number of HFRS cases.
Ordinary least square method (OLS) and geographically weighted regression (GWR): Feature correlation analysis.
**Correlation is significant at the 0.01 level (2-tailed).
*Correlation is significant at the 0.05 level (2-tailed).
The estimation results from the OLS and GWR methods are different because the GWR model includes spatial characteristics in the general linear regression. Compared with the OLS, the GWR model is a nonstationary regression method; OLS has superior performance in the analysis of data with spatial and temporal characteristics, like HFRS trends. 33
Reports in the literature indicate that HFRS in Hubei Province is primarily caused by SEOV, transmitted by Rattus norvegicus. 18 In the present work, Yicheng and Zaoyang were found to be persistent hot spot counties. These counties are mainly covered by farmland, with resident rodent populations. Zuogui, Changyang, and Fangxian were found to be intensifying cold spot counties. Their low coverage of farmland may contribute to the low incidence of HFRS in recent years. Zhongxiang was found to be a diminishing hot spot; its farmland coverage is high in the west and low in the east, so this distribution characteristic has led to a decreasing trend of HFRS incidence in Zhongxiang. Results of local outlier analysis showed that the clustering characteristics in the west of Hubei are mainly in the form of low-low clustering and mainly in the form of multiple distribution clustering in the east of the province; this is in accordance with the increasing proportion of farmland coverage from west to east, as shown in Figure 7. When combined with the results of correlation analysis, the proportion of farmland coverage appears to be positively associated with the emergence and spread of HFRS epidemics in Hubei Province.

Visualization of farmland and water distribution in Hubei Province Blue polygons represent the distribution of water areas; light orange colored areas represent the distribution of farmland in Hubei Province.
The GWR method was used to make predictions based on the spatial and temporal data from 2011 to 2015. In a comparison of Figures 8 and 9, the predicted results for hot spot counties (Zhongxiang, Yicheng) are relatively consistent with the observations, indicating that the selected possible influencing factors are possible determinants of the HFRS epidemics in Hubei Province.

Prediction results for HFRS in Hubei Province, 2016 The color gradient represents the annual predicted incidence of HFRS in thematic maps.

Observed results for HFRS in Hubei Province, 2016 The color gradient represents the annual observed incidence of HFRS in thematic maps.
From the west to the east of Hubei Province, the proportion of water areas increases significantly, which has also led to increasing incidence of HFRS. These results also correspond with the results of correlation analysis. In addition, Peng et al. 39 found that HFRS cases in Anhui Province between 1983 and 1995 were distributed along large water systems (wetlands), such as the Yangtze and Huai river basins. The surrounding water system is beneficial for crop cultivation, which in turn helps to establish thriving rodent populations that are hantavirus hosts. From this study, we infer that the increase of HFRS incidence in Hubei Province from 2011 might be related to the increased population of infected rodents and the surrounding water system. However, in our study, space–time hot spot analysis indicated that counties with sufficient water resources have not always behaved like persistent hot spots. As an example, the annual number of cases in Wuhan from 2011 to 2015 was 33, 12, 0, 5, and 4, respectively. With the combination of disease prevention measures and economic development, HFRS epidemics are effectively controlled in Wuhan County, thereby removing this county from the hot spot clusters.
Human migration is regarded as an influencing factor for the risk of HFRS spread. 17 The space gravity center method has been used to calculate the migration of the space center of HFRS cases and the population. It has been shown that the space gravity center of HFRS cases has consistently changed within a narrow range, as shown in Figure 10. In general, the space gravity center of HFRS cases in Hubei Province has been consistent with that of the overall population from 2011 to 2015. For each year, the change in the range and direction of the space gravity center of HFRS cases and the overall population are also identical. Li et al. found that an increase in human population density could facilitate virus transmission (Xiao et al. 2011). As discussed above, the human population density factor is positively associated with outbreaks of HFRS. In Hubei Province, both human population density and population migration are factors that are significantly related to HFRS outbreaks.

Space weight center moving map and population space weight center moving map of HFRS cases in Hubei Province (2011–2015) Blue lines with arrows connecting the points represent the moving direction of HFRS cases of the space gravity center. Red lines with arrows connecting the points represent the moving direction of the population space gravity center.
Previous studies have reported that the spatiotemporal dynamics of HFRS can be affected by a number of environmental factors, such as moisture (e.g., precipitation, humidity, rainfall) and temperature. 41 In this research, the average humidity and rainfall were found to be positively associated with outbreaks of HFRS. In other lines of research, Xiao et al. 42 demonstrated that rodents thrive in moist environments. Therefore, taken together with the correlation results, it can be inferred that moist climates help to establish a thriving rodent population, which may in turn contribute to HFRS outbreaks.
We recommended initiating a comprehensive prevention and control strategy, a “monitoring, health education, rodent control, and immunization” program, in hot spot areas of HFRS epidemics. Early warning mechanisms are needed in these areas, with continuous routine monitoring of human epidemics, host animal species and density, as well as favorable conditions for virus propagation, especially 1 or 2 months before the peak period of disease outbreaks. We also recommend conducting propaganda campaigns on the prevention and control of HFRS to improve health knowledge among responsible individuals in key areas, and to improve disease prevention awareness and behavior.
Conclusions
With application of emerging hot spot analysis, correlation analysis, and the space gravity center method, the results from our study showed that an HFRS epidemic is more likely to spread in areas with greater farmland or water areas. Population density, human migration, average humidity, and rainfall are also factors related to the incidence of HFRS. These results are in accordance with those obtained in spatiotemporal clustering analysis.
17
The STC method used in this study offered the following advantages:
Based on an STC analysis of disease, the spatial–temporal trend and characteristics of HFRS in Hubei Province could be analyzed from macroscopic and microscopic aspects. Each bin in the STC can be used to analyze and predict every spatiotemporal granularity unit within the space–time range of the entire study period. Based on emerging hot spot analysis and local spatiotemporal outlier analysis, the changes in each spatial unit could be obtained. Combined with the time series of each spatial unit, the degree of potential risk and the development trend of disease could be determined, enabling us to provide constructive suggestions. According to the space gravity center model, the migration characteristics of possible influencing factors could be revealed. Based on the changes and migration of the space gravity center, the distribution and development trend of HFRS and the influencing factors could be analyzed.
In summary, the research method proposed in this paper can be used for the analysis of factors with temporal and spatial distribution characteristics.
Footnotes
Author contributions
Liang Ge and Youlin Zhao conceived and designed the experiments; Liang Ge, Youlin Zhao, Junwei Liu, Honghui Liu, Yijun Zhou, and Ning Wang performed the experiments; Liang Ge, Youlin Zhao, and Xu Ding analyzed the data; Lei Yu contributed reagents, materials, and analysis tools; Liang Ge, Youlin Zhao, Ning Wang, Xu Ding, and Lei Yu wrote the paper; Youlin Zhao, Liang Ge, Junwei Liu, and Honghui Liu revised the paper.
Declaration of conflicting interest
The authors declare that there is no conflict of interest.
Ethics and consent statement
Restricted data are available upon request from the following sources: HFRS data: Hubei Provincial Centre for Disease Control and Prevention. Honghui Liu. E-mail: liuhonghui1987@126.com. Ethical approval was unnecessary in this study. There was no inclusion of nonroutine procedures in this research. Informed consent was also not necessary.
Funding
This work was supported by the National Natural Science Foundation of China, No. 71503068; and National Philosophy and Social Science Fund Key Grant, No.16ZDA046.
