{'fraction_non_alphanumeric': 0.03763350534021361, 'fraction_numerical': 0.016320652826113044, 'mean_word_length': 4.44511238891793, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analysis of the Spatio-temporal Dynamics of COVID-19 in Massachusetts via Spectral Graph Wavelet Theory

The rapid spread of COVID-19 disease has had a significant impact on the world. In this paper, we study COVID-19 data interpretation and visualization using open-data sources for 351 cities and towns in Massachusetts from December 6, 2020 to September 25, 2021. Because cities are embedded in rather complex transportation networks, we construct the spatio-temporal dynamic graph model, in which the graph attention neural network is utilized as a deep learning method to learn the pandemic transition probability among major cities in Massachusetts. Using the spectral graph wavelet transform (SGWT), we process the COVID-19 data on the dynamic graph, which enables us to design effective tools to analyze and detect spatio-temporal patterns in the pandemic spreading. We design a new node classification method, which effectively identifies the anomaly cities based on spectral graph wavelet coefficients. It can assist administrations or public health organizations in monitoring the spread of the pandemic and developing preventive measures. Unlike most work focusing on the evolution of confirmed cases over time, we focus on the spatio-temporal patterns of pandemic evolution among cities. Through the data analysis and visualization, a better understanding of the epidemiological development at the city level is obtained and can be helpful with city-specific surveillance.


I. INTRODUCTION
W HEN writing this paper, the COVID-19 pandemic is still ongoing and has resulted in more than 260 million people diagnosed and more than 5 million deaths worldwide. The COVID-19 spreading was recognized by the World Health Organization (WHO) as a pandemic on March 11, 2020 [1]. As of December 1, 2021, confirmed cases in the United States surpassed 4,810,000 and Massachusetts surpassed 919,000 [2]. With little information on similar past pandemics, collecting mobility, safety, and behavior data related to COVID-19 and learning from the collected data become the key for decisionmakers. The final scale of the disaster had not yet been determined. To make matters worse, the mutation of the virus has led to waves of pandemics, and humankind still faces significant challenges.
Traditional mathematical models use compartmental models to study the transmission dynamics of COVID-19, such as SIR and SEIR models and their variants, see Gao et al. [3], Church [4], Neves and Guerrero [5], Ng and Gui [6]. Miranda et al. [7] construct a hybrid ODE-network model for the COVID-19 pandemic accounting for certain spatial aspects. Some stateof-the-art technologies, such as machine learning, and deep learning, have been introduced in the research of COVID-19. For example, Tang et al. [8] study the interplay of demographic variables and social distancing scores in the deep prediction of COVID-19 cases in the United States. Melin et al. [9] use a self-organizing mapping neural network to cluster countries with similar pandemics so that similar strategies can be used to deal with the spread of the virus. Tat Dat et al. [10] apply wavelet theory and machine learning method to study the evolution of the pandemic in France, Germany, Italy, the Czech Republic, and the US federal states of New York and Florida. Graphs and networks are used to model many realworld problems due to their flexible structure. Li and Mateos [11] conduct a graph Fourier frequency analysis to investigate the county-level contagion patterns of COVID-19. Gao et al. [12] utilize a new graph neural network -STAN to predict both state-level and county-level future number of infected cases. Kapoor et al. [13] examine COVID-19 forecasting using spatio-temporal graph neural networks. However, there is not much work analyzing the spatio-temporal pandemic spread patterns at the city levels.
In this paper, we focus on analyzing the local pandemic spread patterns in the State of Massachusetts (MA), USA, which covers a relatively small area of 21,000 square kilometers. Cars are the primary transportation tools among cities and towns in the state of MA. Therefore, the local spread of the pandemic is mainly affected by the transportation graph structure of the state. Unlike the state-level or countylevel graphs mentioned above, our city-level graph is more suitable for the local epidemiological analysis in MA. In this Massachusetts Route Graph G M R = (V, E), hand-crafted from US Routes, we map cities to vertices on the graph, such that each vertex τ i ∈ V is associated with a set of time-dependent features extracted from biweekly COVID-19 data. The edge between two vertices represents that there is a arXiv:2208.01749v1 [cs.SI] 28 Jul 2022 US Route connecting the corresponding cities. The pandemic spreading between two cities may be affected not only by population movements such as transportation, commuting to work, and tourism; but also by hidden complex factors such as the circulation of COVID-19 virus-contaminated goods on the packaging, etc. Therefore, the traditional method of calculating the edge weight using distance and populations cannot accurately represent the pandemic spread among cities.
This paper further utilizes a deep learning method via graph neural networks -Graph Attention Network (GAT) -to learn the Markov transition matrix P , in order to capture geographical proximity and data features similarity in the transition probabilities. The resulting graph G M AR = (V, E, P ) is called the Massachusetts Attention Route (MAR) graph. GAT has been proved to be efficient for learning the edge weights by adaptive assigning different importance to different nodes through the learning of downstream tasks and data [14]- [19]. We also ranked the top influential cities for pandemic spread based on the attention coefficients learned from the GAT. Moreover, we construct the spatio-temporal dynamic graph model connected by the strong product [20].
To analyze the pandemic signals on the spatio-temporal dynamic graph model, we use a powerful tool -the spectral graph wavelet transform (SGWT). Indeed classical wavelet transformation has played an important role in multiresolution analysis, especially for studying signals containing discontinuities and sharp spikes. However, the construction of wavelets on graphs is rather difficult. Hammond et al. [21] and Shuman et al. [22] construct SGWT, and successfully apply it in studying graph signals and identifying anomalies. The SGWT coefficients contain rich information about the graph signal; however, it remains a big challenge to interpret them properly for non-experts. In this paper, we use a visualization methodology that relies on SGWT and the structure of our graph to enable the visual multiresolution analysis of time-varying COVID-19 data in Massachusetts. We analyze the dynamic patterns of pandemic signals for cities and towns, and provide visualization diagrams based on our spectral graph wavelet analysis results. Moreover, we define an average metric function based on wavelet coefficients to give cities a spatio-temporal ranking synthetically. The top five ranked cities include Springfield, Amherst, Great Barrington, Holyoke, and Sandisfield, all in line with evidence of repeated poor responses to the pandemic, see [23]- [25]. For example, Amherst has reported pandemic peaks during the starting period of semesters of University of Massachusetts (UMass) Amherst, while Holyoke has reported pandemic breakout in its nursing homes. These all make the pandemic situation more severe than those in the surrounding areas. On the other hand, according to our ranking, we identify the five best-performance cities -New Marlborough, Harvard, West Tisbury, Tolland, and Leverett. These are in line with evidence showing their excellent pandemic prevention and control work, [26]- [29], which shows that the search for ambassadors at New Marlborough and the offer of free safety training courses in Harvard, make their pandemic situations much better than that of surrounding cities.
Experts often say that determining which cities and states have had the best response to COVID-19 thus far is tough and unfair [30]. In this paper we hope our study provides a useful evaluation method. The basic framework and main results of this paper are shown in Fig. 1.
In summary, the main contributions of this paper include: • Massachusetts Route graph is constructed by connecting cities and towns based on US Routes or interstate highways. • GAT learns the weights of the MAR graph, which helps calculate each node's vulnerability affected by the COVID-19 infection rate of surrounding nodes. • SGWT is used to analyze the spatio-temporal dynamic model constructed by the strong product. • We design a new node classification method, which effectively identifies the anomaly cities based on spectral graph wavelet coefficients.
The outline of this paper is as follows: In section II, we describe the COVID-19 data of cities and towns in Massachusetts. In section III, we construct the Massachusetts Route graph and learn the pandemic transition probability Matrix by graph attention neural network, then build the spatio-temporal dynamic model of confirmed cases in Massachusetts. Section IV introduces the SGWT and develops the node classification method based on SGWT. In section V, we visually analyze the overall trend of the pandemic by SGWT, make a refined node classification in anomaly cities, and evaluate the city's pandemic prevention work. In addition, we identify the top influential pandemic spreader in the MAR graph. In section VI, we close with a conclusion section.  The dataset available in this paper is composed of population data {N i , i = 1, · · · , 351} and the "two-week-case" time series data {γ i (t), t = 1, · · · , 41} at the city/town level in Massachusetts, from December 6, 2020 to September 25, 2021, for 41 weeks total, which is collected from the official website [31]. We define the biweekly confirmed COVID-19 cases per thousand population (abbreviated as confirmed cases), denoted by the multivariate time series with N = 351 and T = 41.
The multivariate time series data x i (t) will first be used to learn the pandemic transition probability p i,j from city i to city j, in order to better capture the spatio-temporal evolution of the local pandemic spread dynamics. After that, their evolution patterns will be analyzed using the spatiotemporal dynamic graph model, and we are able to identify the most influential spreaders (cities), and detect anomaly cities with unusual pandemic spread patterns.

SPATIO-TEMPORAL DYNAMIC MODEL
In this section, we first construct the Massachusetts Route graph G M R = (V, E) based on US Routes, then learn the attention graph weight P from the COVID-19 data using GAT to get the Massachusetts Attention Route (MAR) Graph G M AR = (V, E, P ). The strong product is utilized to construct a spatio-temporal dynamic graph model.

A. Construction of the Massachusetts Route graph
We take the longitude and latitude coordinates of the cities and towns as the vertices τ i ∈ V from the website [32]. The detailed list of cities and their corresponding vertex ID is included in the Appendix. An edge e i,j ∈ E represents the existing US Routes or interstate highways go through the administrative regions of the cities τ i , τ j ∈ V successively. We also connect the cities of Barnstable (No.21), Nantucket (No.197), and Oak Bluffs (No.221) because of the heavy sea transportation among them. The number of cities and towns in our study is |V | = N = 338 1 . Let A ∈ R N ×N be the adjacency matrix of the graph G M R , such that a i,j = 1 if e i,j ∈ E, and a i,j = 0 otherwise. Each vertex has signals x i ∈ R T , γ i ∈ R T and N i ∈ R defined as in formula (1).
The size of vertices in Fig. 3 is proportional to the city population. For example, Vertex No. 36 represents metropolitan Boston, with the largest population in Massachusetts. Since the cities within the range between 71 • E to 71.5 • E and 42 • N to 42.8 • N are dense, to improve the visualization in Fig. 3, we use the downsampling method [33] by hiding some nodes in the Massachusetts Route graph. Specifically, we hide the nodes with negative signs in the nodal domain of the largest eigenvector of the Laplacian of the Massachusetts Route graph G M R .

B. Massachusetts Attention Route graph learned by graph neural network
The pandemic spreading between two cities may be affected by many complex factors, including transportation, tourism, commuting, and other hidden factors such as the circulation of 1 There are 13 small towns located in rugged areas without any highway or heavy sea transportation. We do not consider them in our following study. COVID-19 virus-contaminated goods. Due to the complexity of the pandemic, the traditional weighting method by distance is not a good choice in modeling the pandemic spread between cities (towns). In order to accurately capture the pandemic transition probabilities between cities, we utilize a deep learning method via GAT, to construct the transition probability matrix P = {p ij }, based on the importance of geographical proximity on the Massachusetts Route graph G M R = (V, E) and data feature similarity. The resulting new directed graph with self-loops is called the Massachusetts Attention Route (MAR) graph, denoted as G M AR = (V, E, P ). The edge set E includes E, as well as self-loops.
GAT is a powerful deep learning method proposed in [14], which has state-of-the-art performance on node classification and link predictions. Instead of a simple mean aggregation, the GAT learns to use weighted summation to aggregate the features of neighboring nodes to update the features of the nodes, and automatically adjust the weight through downstream tasks [14]- [18]. Given N multivariate time series x i (t), for i = 1, · · · , N , t = 1, · · · , T . Then x i ∈ R T can be viewed as the T -dimensional input features of the node τ i . We use edge classification as the downstream task, for our GAT neural network, and use the normalized attention matrix in the last layer of the neural network as our Markov transition To improve the model capacity and stabilize the learning process, GAT uses multi-head attention. Consider a multi-head be the k-th attention weight vector. The normalized attention coefficients α k i,j capture the importance of τ j 's features to τ i 's, expressed as where N i denotes the set of first-order neighbours of node τ i (including τ i ) and denotes the concatenation operation, and the activation function is LeakyReLU 2 . Let A K-head GAT layer aggregates nodes features across neighborhoods by where ELU is an exponential linear unit 3 .
To learn the attention matrix, we use edge classification [15] as a downstream task. Let be the positive sample set, where a i,j is the element of the adjacency matrix A of the graph G M R . For any integer l ≥ 1, be the negative candidates set. The negative sample set M − is randomly generated from the negative instances in S − with a similar size to M . The GAT neural network is trained by minimizing the following cost function: Here q ij is the probability of whether there is an edge between τ i and τ j , where denotes Hadamard product, x i is the output feature of the last layer of GAT. θ is a learned parameter vector and the dimension is the same as x i . In our setup, there are two GAT layers. Each layer takes an appropriate output dimensionality to improve the performance of the model [34]. The first layer is a K-head attention layer with K = 7, and the input feature is the graph signal i.e. we set O k = 122 features each (for a total of 854 features). The second layer is a single-head attention layer with the 88dimensional output features x i . During the training phase, we divide the sample set M ∪ M − into training set, validation set and test set according to the ratio of 6 : 2 : 2. We also make sure to keep the same number of positive and negative samples in the test set. The learning rate is set as lr = 0.005. The early stopping strategy is applied to the validation set to avoid overfitting, with the patience set to 100 epochs. The accuracy of edge classification of G M R without self-loops in our result is 0.9499.
The transition probability p ij is defined as the GAT coefficients of the last GAT layer, i.e.
, whereâ and W defined above are the learned attention weight vector and the learned parameter matrix of this singlehead attention layer. One advantage is that the self-attention coefficient p i,i > 0, which also explains the self-influence of the pandemic in a city. Consequently, the obtained MAR graph is a graph with self-loops. Fig. 4 shows the basic framework of this method.

C. Spatio-temporal dynamic graph model
To perform the spatio-temporal node classification using the time-dependent data {x i (t), i = 1, · · · , N, t = 1, · · · , T } defined in (1), we introduce the concept of graph product. The three well-known graph products are the graph Cartesian product, graph Kronecker product, and graph strong product, see Fig. 5. These have been powerful in studying time-vary data on graphs, see [35]- [38].  The graph Cartesian product only adds extra edges between nodes on the adjacent graph time slices in series, which often weakens the influence of spatial neighbors in the study of SGWT for graph signals. Here we adapt the graph strong product, which also adds temporal edges between pairs of onehop neighbors in its adjacent graph time slice.
Let H be a path graph consisting of T nodes. The strong product graph can be defined by the product of the Massachusetts Route graph G M R and the path graph H. We assign weights to the strong product graph based on the transition probability matrix P of G M AR . The resulting directed graph G without self-loops is the spatio-temporal dynamic model. More precisely, the graph time slices {G t = (V t , E t , P t ), t = 1, · · · , T } are copies of graph G M R by assigning weights P t = P − diag(p 11 , p 22 , · · · , p N N ). The direction of the temporal edge (black edge in Fig. 6) is from G t to G t+1 , which represents the time-evolution of the pandemic among these cities. The weights of the new temporal edges are carried from P . Now, the strong product graph can be represented as G = (V, E, W) with V = {τ i,t , i = 1, · · · , N ; t = 1, · · · , T }, and E = (E 1 ∪· · ·∪E T )∪(E 1,2 ∪· · ·∪E T −1,T ), where E t,t+1 is the collection of temporal edges connecting G t and G t+1 . The new weight matrix W is not normalized anymore, as we have added extra temporal edges, which carry the same weight as the spatial edges. More precisely, the edge (τ i,t , τ j,t+1 ) has the same weight p ji as the spatial edge (τ i,t , τ j,t ), see Fig.   6. Note that W is rather sparse, as we only added temporal edges between two adjacent graph time slices. Fig. 6. Construction of spatio-temporal dynamic graph based on MAR graph. The black edge in the spatio-temporal dynamic graph represents the temporal edge, and the blue edge is the spatial edge inherited from the MAR graph without the self-loop.

IV. NODE CLASSIFICATIONS ON THE SPATIO-TEMPORAL DYNAMIC GRAPH G A. Brief introduction to the Spectral Graph Wavelet Transform
Let G = (V, E, W) be the spatio-temporal graph constructed in the above section. The graph signal X on V with X(τ i,t ) = x i (t) defined on (1) represents the confirmed cases of city i in the t-th week. Clearly, the cardinality of |V| =N := N T . The graph Laplacian is defined by L = D − W, where the degree matrix D = (d ij ) is a diagonal matrix with entries d ii = k w ik . The graph Laplacian L is not symmetric as our pandemic graph is directed. We construct a symmetric matrix by L = (L+L * )/2, which is the averaging of the graph Laplacian and its transpose. L hasN non-negative, real-valued eigenvalues Λ = {λ 1 , λ 2 , · · · , λN } with λ 1 = 0, and the corresponding normalized eigenvectors {u l , l = 1, · · · ,N } form an orthonormal basis for the Hilbert space. The graph Fourier transformation of X can be written aŝ A graph spectral filter, or kernelĝ : Λ → R is a function defined in the spectral domain λ l ∈ Λ. A dictionary {ĝ m } m=1,2,··· ,M is a set of graph spectral filters constructed to detect different frequencies of signals.
In this paper, we use the dictionary constructed by Hammond et al. [21]: where the wavelet kernel is a bandpass filter defined on the Fourier domain, the stretching scale s 1 , s 2 , · · · , s M −1 are logarithmically sampled between s 1 = 1/λN and s M −1 = 40/λN , the scaling function is a low-pass filter detecting the low frequency signal with the parameter b = max λĝ (λ). Such dictionary is designed to evenly cover the graph spectrum domain. Given a graph signal X and the dictionary {ĝ m } m=1,2,··· ,M , the spectral graph wavelet coefficient However, the above SGWT is rather expensive for large graphs, as the resulting computational complexity is of order O(N 3 ). In this paper, we apply the fast spectral graph wavelet proposed in [21] to overcome this difficulty. The fast spectral graph wavelet transform based on Chebyshev polynomials approximation [39] shows that the graph wavelet coefficients in (4) can be approximated by where K m is the number of truncating terms, and c m,k = 2 π π 0 cos(kθ)ĝ m (cos(θ))dθ.
Moreover, T k is the shifted Chebyshev polynomials with the domain of [0, λN ], such that T k (x) satisfies the recursive formula The computational cost to approximate the wavelet coefficients is order O(K m |E| + K m |V|).
In this paper, we take K m = 40 in (5). It takes 0.5 seconds on a laptop with a 2.3 GHz 8-core Intel Core i9 processor. Without Chebyshev approximation, it costs 650 seconds.

B. Node Classification Using SGWT coefficients
Given our spatio-temporal dynamic graph G = (V, E, W) together with the COVID-19 confirmed cases X, one question is how to provide a high-quality classification of nodes according to the pandemic spreading patterns? We refer to this as the "node labeling problem" with the understanding that the node classification problem can be abstracted as providing a classification label for the graph structure. Here we will use the information already encoded in the SGWT coefficients to help us predict labels.
Let T = 41 be the number of time slices in the spatiotemporal dynamic model from December 6, 2020, to September 25, 2021. Let M = 8 be the number of graph spectral filters. Using the SGWT on the graph signal X, one gets â N × M wavelet coefficient matrix with entires (5). We then define node classification criteria based on the multiresolution information carried by SGWT coefficients. We will classify the nodes by the torque value obtained by some weighted average of the SGWT coefficients.
Step 1: Data processing of wavelet coefficients.
Firstly, we use a data scaling method to ensure for each wavelet coefficient that guarantees that they are on the same scale. This method was introduced by RobustScaler [40] in the machine learning library Scikit-learn as follows: where R(m) is the interquartile range of |W X (m, τ )| with respect to τ . The interquartile range is the difference between the first and third quartile of a data set. The quartiles make the RobustScaler ignore data points that are very different from the rest (like machine errors).
Logarithmic normalization does not change the order of the values, but it can significantly reduce the effect of maximum abnormal value [41].
Step 2: Spatial node classifications based on torque values. In order to use the SGWT coefficients to classify the nodes, we need to define the torque value function ϕ : V → R by: High torque values imply a high frequency of signal changes [38]. More precisely, if a vertex τ (i, t) has a relatively high torque value, then there are anomalies in the pandemic spread patterns of the i-th city compared to its surrounding cities at the time t. Let d = max τ ϕ(τ ) − min τ ϕ(τ ). By splitting the interval [ϕ min , ϕ max ] into 5 equal intervals, the vertex space V can be divided into five disjoint sets {V 1 , V 2 , V 3 , V 4 , V 5 }, which represent the spatio-temporal changes of the x i (t) (in terms of the SGWT coefficients) as low frequency, mid-low frequency, uncertain, mid-high frequency, and high frequency, respectively. The detailed node classifications are defined in Table I. The spatial node classification method is demonstrated in Fig. 7, for the 2-week period of Jan.31-Feb. 13, 2021. An efficient node classification should be able to capture the abrupt changes of the graph signal, both spatially and temporally. To make the illustrations more explicit, we use our node classification method to investigate the city of Amherst during the 13 weeks from December 16, 2020, to March 13, 2021. Fig. 8 shows that the confirmed cases in Amherst were significantly higher than its neighbors from January 24, 2021 to February 27, 2021. Fig. 9 demonstrates the different temporal stages of the node classification value of Amherst during the 13 weeks. Compared with Fig. 8, we can see that when the City of Amherst has similar pandemic patterns as the surrounding cities, the wavelet coefficient represents low frequency (Dec. 6-Dec. 19). On the other hand, as the confirmed cases in the surrounding cities start to rise, the wavelet coefficient of Amherst represents high frequency (Dec. 20-Jan. 2), which means that the pandemic cases of Amherst have different patterns than its neighbors. Our node classification method accurately captures the abnormal pandemic spread patterns of Amherst during the period from January 24, 2021, to March 6, 2021, see Fig. 8 and Fig. 9. In this period, the University of Massachusetts Amherst partially opened its campus, and many students returned to Amherst [42].
This method can effectively identify the pandemic spread anomaly patterns of any specific period, which should be valuable for further research.

V. VISUALIZATION ANALYSIS OF COVID-19 PANDEMIC SPREAD PATTERNS
In this section, we focus on the visualization analysis of the COVID-19 pandemic spread patterns on the spatio-temporal dynamic graph model. We first conduct an overall analysis of the COVID-19 data in Massachusetts. Then we identify the supper pandemic spreaders (cities). Finally, a refined node classification is performed in V 4 and V 5 , where the pandemic spread patterns have abrupt spatio-temporal changes.

A. Overall graph classifications
To have an overall temporal analysis of the pandemic evolution patterns for cities in Massachusetts, we define a graph classification method for the spatio-temporal dynamic graph G, which enables us to classify each graph time slice for {G t , t = 1, · · · , T }.
First, we need to calculate the node classification distribution in G t , which can be obtained by counting the frequencies of nodes in G t that occur in each of the five classes V j , j = 1, · · · , 5: where I Vj is the indicator function of the set V j . The distribution {σ j t , j = 1, · · · , 5} can not be used to directly classify the graph G t , as one has to weigh in the distributions in other time slices to well capture the role of each frequency class. We adopt the method used by Dal Col. etc. [38], by first taking σ j max = max{σ j 1 , σ j 2 , · · · , σ j T }. Then we define the class of G t as V rt , with r t = arg max j∈{1,2,3,4,5} (σ j t /σ j max ).
Next we conduct a temporal graph classification to get an overall visual analysis of the spatio-temporal pandemic spread patterns in Massachusetts, which is demonstrated in Fig. 10. As we can see from Fig. 10, after Thanksgiving, 2021, Massachusetts saw a substantial spike in confirmed cases as there were more group gatherings. The confirmed cases in   (7), and the spatial node classification distribution σ j i given by (6) (bottom).
Massachusetts reached their peak during the Christmas holiday period. We can infer that multi-site outbreak patterns exist because the temporal and spatial changes of confirmed cases in many cities are different from those in the surrounding cities. The cases declined sharply after mid-January due to restricting gatherings, maintaining social distance, and advocating masks, then stabilized in late February (mid-low frequency). Massachusetts entered the fourth phase of the state's reopening plan on March 22, 2021. With the relaxation of pandemic restrictions, there was a small-scale rebound of confirmed cases with the multi-point pattern at the end of March (midhigh frequency). The numbers began to decrease from April to July due to the increased vaccination rate.

B. Spatial ranking for super-spreader cities
Identifying influential spreaders (cities) in a pandemic network is a core question to prevent pandemic spreading. Many different centrality measures tied to the network topology have been introduced to find the central nodes. However, the importance of node features is mostly ignored in these works. The centrality of a node in a pandemic network should depend on two features: its local influence on the nodes in its one-hop neighbor and its global influence on the nodes belonging to higher-order hop neighbors. Taking advantage of the attention coefficients learned from the GAT, we introduce a new score function that allows us to redefine the classical centrality measure by considering their importance to the neighbors' features.
Considering the m-hop neighbor of node τ i in the MAR In our practice, we define the influential score of city i by as the total capability of city i to spread the virus to other cities.
The spatial ranking of cities according to their influential scores is shown in Fig. 11. One can see that the top ten cities with the largest pandemic spread (influential) scores (in dark red) are Boston, Worcester, Springfield, Otis, Fitchburg, Newton, Halifax, Dedham, Montague and Palmer. It's not surprising that Boston, Worcester, and Springfield are in the top three, given these cities' large populations and high mobility. On the other hand, we are able to identify small cities, which are otherwise difficult to identify without using our advanced methods. For example, the City of Otis has a small population but ranked fourth on the list. It is remarkable that the City of Wayland has a larger population and the same number of connected US Routes compared to Otis, but it ranks much lower. This has been confirmed by news in [43], as well as the fact that Otis has a ski resort. Using this new ranking metric, we can correctly identify cities at risk and those more volatile to pandemics and thus need to invest more social resources to better control the pandemic.

C. Refined node classification in anomaly cities
Given the COVID-19 data in Massachusetts, another important question is to identify cities with anomaly spread patterns during the pandemic period. More precisely, this is equivalent to finding patterns in data (i.e., attribute values or changes in the values over time) that are significantly different from that of the spatio-temporal neighbors. These non-conforming patterns are often referred to as anomalies. The typical anomaly detection method concentrates on time series, or static graph signals, which could not be applied to our case. We aim to design a new method that enables us to detect anomalies in COVID-19 time-series data belonging to multiple entities (cities) and explain the anomalies to domain experts in a comprehensible manner.
SGWT is a powerful tool for finding anomalies for spatiotemporal graph signals. For example, the set V 5 can be counted as one anomaly class, both spatially and temporally. However, in V 5 , by examining the distribution of the SGWT coefficients more carefully, we find that there are situations in which there is a node with a higher or lower signal value than that of the surroundings. As shown in Fig. 12, there are two dramatically types of vertices in V 5 , as indicated above.
Below we say the nodes in V 4 ∪ V 5 are anomaly nodes. The corresponding city is the anomaly city. We use a refined node classification method to identify these two abnormal situations in anomaly cities. Let N We first define an anomaly metric: which is the relative ratio of the data information of city i at time t compared with the averaging data obtained from its one-hop neighbors. Then we define the anomaly score (a-score) a(i, t) to the nodes in the set V 4 ∪ V 5 according to Fig. 13. According to the a-score, we divide the V 4 ∪ V 5 into five groups as in Fig.  13. A city with high a-scores presents at risk or volatile in the pandemic. It should ring an alarm for local city administrations, as the city has the potential to increase infection spread to its neighboring cities. Earlier identifying these cities can prevent them from becoming the outbreak cities of a pandemic. On the other hand, identifying cities with low a-scores can be valuable in controlling virus spread, as they provide good models for other cities.

1.31-2.13
Fig. 14. Refined node classification for the 9th week, from Jan.  Fig. 14 depicts the data visualization plots for cities with anomaly patterns from January 31 to February 13. We can quickly detect Amherst with a high a-score, which is also confirmed as one of the most significant COVID-outbreak cities in Fig. 9.

D. Spatio-temporal ranking of super-spreader cities
The a-score defined in Fig. 13 depends on both spatial and temporal coordinates. To have an overall picture of the pandemic evolution behavior of these cities and evaluate the urban pandemic prevention and control, we propose a data visualization plot to explain the level of anomalies by the anomaly score. Below, we use a-scores in Fig. 13, to provide a new city ranking, by taking averages of its a-scores over the 41-week period for each city:  One can see that cities with higherā(i) values are more likely to become super-spreader (cities) for any future pandemic. A by-product of using theā(i) value is that, one could evaluate the city's pandemic prevention and control capabilities. In Fig.  15, the darker the red, the cities have the less successful pandemic prevention and control. On the other hand, the darker the blue, the cities are more successful in pandemic prevention and control. The top five least successful cities are identified as Springfield, Amherst, Great Barrington, Holyoke, and Sandisfield; as well as the top five most successful cities are New Marlborough, Harvard, West Tisbury, Tolland and Leverett.
The detailed visualization analysis for these cities is included in Fig. 16. This visual view can quickly locate the outbreak time and duration of abnormal phenomenons. To make effective comparisons, we also plot the pandemic time evolution data of Springfield, Amherst, New Marlborough and Harvard, see Fig. 17 -Fig. 20.
Our analysis identifies Springfield as the least successful city using data from December 6, 2020, to September 25, 2021. Fig. 16 shows that most of the confirmed cases in Springfield surpass the values of the surrounding areas during the entire 41-week period, which can also be confirmed in Fig.  1  17. A recent report shows that Springfield is a center for the Delta outbreak [23]. This confirms that the city has a much higher risk of spreading new pandemic variants to other cities and may need more help from administrations to control the pandemic. Amherst is ranked as the second least successful city for 41 weeks from December 6, 2020 -September 25, 2021. Fig.  16 indicates that Amherst has much higher confirmed cases than neighbors from January 24 to March 27, 2021, which is consistent with Fig. 18. Indeed, we show that Amherst ranks as the first pandemic spreader during a two-week short periodic from January 31 -February 13, 2021, see Fig. 7. By digging into the news of this period, we find that the University of Massachusetts Amherst started the spring semester on February 1. Because many students moved into dorms or apartments in Amherst, there had been more than 398 active cases during that period. Contact tracing data record [44] has shown that some students failed to follow social distancing and mask protocols in social and residential settings, promoting pandemic virus transmission. Subsequently, the university further strengthened its management standards. For example, students who do not wear masks will not be allowed on campus, and large group gatherings were not encouraged. The pandemic spread was effectively brought under control at the end of February 2021.
At the beginning of the fall semester (the end of August 2021), a similar abnormal situation appeared again in Amherst, see [24], [45]. New Marlborough and Harvard are the top two most successful cities for the 41 weeks from December 6, 2020, to September 25, 2021. Fig. 16 illustrates that the two cities have much lower confirmed cases than those in the surrounding cities most of the time. Fig. 19 and Fig. 20 describe the time Fig. 19. The 41-week evolution of COVID-19 confirmed cases for New Marlborough and its neighboring cities from Dec. 19, 2020 to Sept. 25, 2021. evolution of the top two most successful cities, which have good and stable spread patterns during the pandemic, and are consistent with the situation illustrated in Fig. 16. It indicates Fig. 20. The 41-week evolution of COVID-19 confirmed cases for Harvard and its neighboring cities from Dec. 19, 2020 to Sept. 25, 2021. that these cities might have adequate prevention and control measures. By digging into the news of this period, the Health Board in New Marlborough looked for a temporary parttime COVID Ambassador to provide COVID-19 prevention guidance to citizens, businesses and community groups [26], and held several sessions for discussions and updates on COVID-19 [46]. Furthermore, New Marlborough has provided free masks to residents since August 20, 2020 [27]. These measures allow New Marlborough to remain unaffected even when its neighbors sharply all had a sharp increase in cases, see Fig. 19. Harvard offered free safety training courses for workers [28]. Citizens volunteered to help the Board of Health open a local vaccine clinic [29]. These measures may be the reason for the success of Harvard's pandemic prevention.
In Fig. 16, we show the five most successful and least successful cities. By changing the time range of the formula (9), we can also find successful and least successful cities for any period. As a by-product, our visualization analysis is based on the refined node classification for our spatio-temporal dynamical graph, so it captures correlations of pandemic evolution patterns among these cities. Take Springfield as an example, see Fig. 16. Because the spatio-temporal outbreak patterns of Springfield and Holyoke are similar, we can infer a significant correlation between them. Indeed, we know that they both have large traffic flows and are close to each other, see Fig. 3, thus the patterns of their pandemic outbreaks should be highly related. Our method can quickly help experts locate cities with similar traits and apply insights from a city to inform prevention and control strategies in similar areas.

VI. CONCLUSION
This paper draws a Massachusetts Route graph based on the latitude and longitude of 351 cities and towns and the main traffic routes passing nearby. The pandemic transition probability is learned through semi-supervised deep learning, based on the Graph Attention Neural Network. We construct a spatio-temporal dynamic model by strong-product for the time-varying data on the graph, which can better capture the spatio-temporal pandemic evolution patterns. With the help of SGWT, the COVID-19 pandemic spread patterns in Massachusetts from December 6, 2020, to September 25, 2021, are analyzed and visualized. In addition to the overall analysis of the temporal and spatial evolution of the pandemic data, we also identify cities with strong pandemic spread influence, such as the town of Otis. We also construct a new anomaly indicator to classify cities with abnormal pandemic situations, indicating cities at risk of virus transmission, and identifying cities that are likely to be infected in the next step. The five cities identified with the least successful pandemic prevention are Springfield, Amherst, Great Barrington, Holyoke, and Sandisfield; the five cities with the most successful pandemic prevention are New Marlborough, Harvard, West Tisbury, Tolland, and Leverett. These obtained results not only offer descriptive insight for strategizing purposes in combating the COVID-19 outbreak in Massachusetts, but also can be used to protect those high-risk cities for the next possible round of pandemics. Further, it also provides a framework for evaluating regional pandemic prevention work in the COVID-19 pandemic. It helps the policymakers identify cities with better prevention in the pandemic and extract insights from the successes.
However, in this paper, we only analyze cities and towns within Massachusetts and do not consider the impact on other states. In our further study, we will take our focus on the influence of neighboring states to improve our method. In addition, when calculating the transition probability matrix using GAT, due to the limited features of our data, we take the time series of confirmed cases every two weeks and obtain the transition probability of the comprehensive pandemic evolution for 41 weeks. Our further research will try to overcome the limited features of the data and train the network on the spatio-temporal dynamic graph to obtain the time-varying transition probability matrix.