{'fraction_non_alphanumeric': 0.036714216940735296, 'fraction_numerical': 0.01385968023918717, 'mean_word_length': 4.734960767218832, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 5, 'lorem ipsum': 0, 'www.': 10, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physically-interpretable classification of biological network dynamics for complex collective motions

Understanding biological network dynamics is a fundamental issue in various scientific and engineering fields. Network theory is capable of revealing the relationship between elements and their propagation; however, for complex collective motions, the network properties often transiently and complexly change. A fundamental question addressed here pertains to the classification of collective motion network based on physically-interpretable dynamical properties. Here we apply a data-driven spectral analysis called graph dynamic mode decomposition, which obtains the dynamical properties for collective motion classification. Using a ballgame as an example, we classified the strategic collective motions in different global behaviours and discovered that, in addition to the physical properties, the contextual node information was critical for classification. Furthermore, we discovered the label-specific stronger spectra in the relationship among the nearest agents, providing physical and semantic interpretations. Our approach contributes to the understanding of principles of biological complex network dynamics from the perspective of nonlinear dynamical systems.

Complex biological systems are modelled as a collection of the discrete elements that non-linearly interact. Dynamic processes for such complex systems are of great interest in a variety of scientific fields, such as sociology 1,2 , epidemiology 3,4 , neuroscience 5,6 and physics 7,8 . These systems are often represented as networked systems using graphs, which have several problems such as in classification 9,10 , prediction 11,12 and control 13,14 . Among various networked systems, dynamic networks of collective motions [15][16][17] , in which nodes and links are agents and their relationships (e.g. based on their positions), pose several challenges. These challenges arise in changes of the relation among agents (i.e. the structure of a graph) in transient and complex ways, with the exception of moderate changes such as diffusion or contagion 12 . In the network dynamics of collective motions, if the governing equation is given such as mathematical models (e.g. cellular automata 18 and coupled dynamical systems 19 ) and controlled systems 20 , these are often modeled as graph dynamical systems (GDSs). However, for many biological collective motions 21,22 , the relation among agents transiently and complexly changes according to different situations, with the exception of the simply modelled collective motions 23,24 . These biological network dynamics can be regarded as the network dynamics or GDSs of collective motions represented by graph sequence data. To address this problem, conventional approaches have basically computed the properties of a graph in each temporal snapshot 5,25 or in a temporal sliding window 26 of the sequence data. However, for the understanding of network dynamics, these approaches are difficult to directly extract the dynamical properties, i.e. physically-interpretable information about the dynamics such as frequencies with decay/growth rate and the corresponding spatial (network) structures. The fundamental question addressed here is how complex collective motion networks should be classified based on the physically-interpretable dynamical properties.
The motivation of this paper is to understand the principles of network dynamics of biological complex collective motions by directly extracting the dynamical properties of the network in a data-driven manner. As a method of describing nonlinear dynamical systems with a global mode by the direct extraction of dynamical properties, operator-theoretic approaches have attracted attention in fields such as applied mathematics, physics and machine learning. One of the approaches is based on the composition operator (often referred to as the

Results
Graph DMD framework. Here, we briefly review Koopman spectral analysis, which is the underlying theory for various DMDs, and then describe the Graph DMD framework. First, we consider a nonlinear dynamical system: = , where x t is the state vector in the state space  ⊂ M p with time index T N ∈ = t : 0 . The Koopman operator, which we denote by K, is a linear operator acting on a scalar observable function M → g: where f g denotes the composition of g with f 27 . That is, it maps g to the new function f g . We assume that K has only discrete spectra. Then, it generally performs an eigenvalue decomposition: is the j-th eigenvalue (called the Koopman eigenvalue) and ϕ j is the corresponding eigenfunction (called the Koopman eigenfunction). We denote the concatenation of scalar functions as  , where ψ j is a set of vector coefficients called the Koopman modes. Through the iterative applications of K, the following equation is obtained: are observed values for a vertex = … i m 1, , . They are given by = y g x : ( ) m m is an adjacency matrix, whose component A i j t , , represents the weight on the edge between V i and V j at each time t. For example, the weight represents the relation (e.g. a function of distance) between moving agents in multi-agent systems.
For the formulation of Graph DMD, we propose a more straightforward formulation with dependent structure among observables than the previous formulation based on vector-valued reproducing kernel Hilbert spaces (RKHSs) 33 . Regarding the details of the analysis, the connection with Graph DMD (such as the relation between A t and g c ) and the reason why it can visualise the relation between elements, see Materials and Methods. For the practical implementation of the spectral decomposition, a modified tensor-based DMD 33,42 , which is a generalised DMD for the application to tensor data, is applied to the adjacency matrix series.
In Graph DMD, a sequence of the matrices A t or an order-3 tensor A t :,:, for τ = … t 0, , is given, where colons are used to indicate all elements. In a similar way of basic DMD procedure (see Materials and Methods), we define tensors X Y , such that = X A by matricising after tensor-train decomposition of Y, where s 2 is a tensor-train rank. Note that this is similar to SVD in the matrix form, but SVD and this matricisation after tensor-train decomposition with maintaining the tensor structure are completely different. After that, in a similar way of basic DMD procedure, we then define a matrix 1 , where M* is the Hermitian transpose of M and N † is the pseudo-inverse of N. Thereafter, we perform eigendecomposition of F and obtain eigenvectors w j and eigenvalues λˆj for = … j p 1, , . The latter is the estimated Koopman eigenvalues called Graph DMD eigenvalues. Finally, we obtain the spatial in the inputted adjacency matrix form, which are called Graph DMD modes. In summary, a sequence of adjacency matrices is decomposed into spatial coefficients and temporal dynamics ( Fig. 1b): Examples of Graph DMD for ballgame data. Next, we show a representative example of Graph DMD results in Fig. 2. Here we used player-tracking data from actual basketball games. Position data was composed of the horizontal Cartesian positions of every player and the ball on the court, recorded at 25 frames per second. We analysed an attack-segment defined as the period that begins when all players enter the attacking side of the court and end before a shot. The input of Graph DMD is adjacency matrix series A t visualised in Fig. 1a and its elements are: where i and j indicate players and a ring (i.e. a goal as geometric information), ||·|| 2 is the Euclidean norm and σ is a coefficient for adjusting the value of A i j t , , , which represents a proximity with contextual meaning in this case (for details, see Materials and Methods). Furthermore, since the order of i is not uniquely determined in general, we sorted the adjacency matrix in order of nearest attackers and defenders from the ball in each time stamp, to Scientific RepoRtS | (2020) 10:3005 | https://doi.org/10.1038/s41598-020-58064-w www.nature.com/scientificreports www.nature.com/scientificreports/ provide the order of player to the semantic information. Here, we denote the first to the fifth nearest attackers and defenders to the ball as A1, ..., A5 ( = … i 1, , 5) and D1, …, D5 ( = … i 6, , 10), respectively (for the ring, = i 11). Also, we denote the relationship between two players, e.g. A1 and D1 as A1-D1.
Firstly, we obtained DMD eigenvalues via Graph DMD. As illustrated in Fig. 2b, all eigenvalues appear to be on the unit circle. Eigenvalue λˆj is transformed into temporal frequency ω j such that ω λ π where Δt is a time interval of the discrete time system (i.e. Δ = t 1/25). Then, for each eigenvalue, we obtained the time dynamics in Fig. 2c-e and Graph DMD modes in Fig. 2f-h. For the sake of clarity, only three pairs of the dynamics and modes showing the smallest frequency in this order are presented. In Fig. 2c-e, Graph DMD extracted the dynamics with approximately one, two and three cycles in a sliding window, indicating that two players with these modes approached and left once, twice and three times, respectively. For the Graph DMD modes, we visualised the strength of the spectrum (i.e. coefficient) for each time dynamics. Note that DMD computes complex-valued DMD modes and time dynamics, but here we only present the absolute value and real parts, respectively. In Supplementary Text 2, we describe the selection of parameters for Graph DMD (e.g. tolerance, temporal window size and cutoff frequency) and quantitatively validated the applicability of Graph DMD to the sport data in terms of the reconstruction error.

Classification in various global collective behaviours.
Next, we indicate that our methods had better classification performance than most of the existing methods in Table 1 and Fig. 3. As validation, we performed two classification tasks with different global collective behaviours: the team-defence (zone or person-to-person defence) and team-offence (offence with and without screen-play) recognition tasks. By definition, zone and person-to-person defence are exclusive (i.e. players guard their area and opponent players, respectively). However, they are actually mixed based on different situations whereas experienced people can distinguish them. A screen-play is the basic and minimal strategic cooperative play in basketball 44 , in which an attacker stands on the course of defence player like a'screen' and prevents the defender from defending another attacker in a legal way. Additional details are provided in Materials and Methods.
After Graph DMD, we performed classification using three different approaches. The first two approaches create feature vectors and the third automatically extracts the features via a neural network approach. The first approach is simply to vectorise the Graph DMD modes (denoted GDMD spectrum) to directly reflect the node (i.e. player) information. For the team-defence recognition task, we used the elements of the modes regarding the relations among defenders, those among attackers and defenders and those among defenders and ring (as geometric information) as a feature vector. For the team-offence recognition task, we used those among attackers and among attackers and defenders (without geometric information). Selection and elimination were based on prior knowledge of the sport. The final two approaches computed graph features using existing methods: the second is graph Laplacian eigenvalues 39 as a baseline feature of a graph 45 (denoted GDMD Laplacian) and the third is deep graph convolutional neural networks 46 (denoted GDMD GCN) as a recent promising method for classifying the graph data. The second approach extracts the graph topology without node information, and the third one retains more node information and learns the graph topology.
Due to the highly-strategic (or contextual) behaviours in this experiment, we hypothesised that, in addition to physical information, the contextual node information (i.e. the nodes themselves in the first approach) is more important for the classification than the graph features (i.e. the second and third approaches). For all classification methods with the exception of the neural network approaches, we adopted logistic regression as a simple linear binary classification model.
As comparable methods for classifying this data, we adopted four methods. The first is a method using vectorised basic DMD modes 47 (denoted DMD spectrum) as a baseline of DMD approaches. The second is the Koopman spectral kernel 40 using DMD with reproducing kernels 41 as the existing method for classifying the collective motion dynamics 38 (denoted KDMD spectral kernel). For these two methods, the input data must be a matrix; thus we did not utilise the structure of the graph sequence data. The third is a hand-crafted feature as a simple baseline method, consisting of the vectorised temporal average, maximum and minimum values of the elements of the input adjacency matrix series. Fourth is an advanced neural network approach for classifying graph sequence data called spatio-temporal GCN 36 as a recent promising end-to-end method. More information on the selection and the details of these methods are provided in Supplementary Text 1.  Table 1. Classification performance of seven methods for two recognition tasks. Accuracy (Acc), area under the curve (AUC) based on receiver operating characteristic (ROC) curve and F-measure are indicated. Overall, classification performances in GDMD spectrum were higher than those in most of remaining methods. (2020) 10:3005 | https://doi.org/10.1038/s41598-020-58064-w www.nature.com/scientificreports www.nature.com/scientificreports/ We investigated classification performance in Table 1 in terms of accuracy, the area under the curve (AUC) based on receiver operating characteristic (ROC) curve in Fig. 3a,b, and F-measure, which is the trade-off between recall and precision (the curve is shown in Fig. 3c,d). Overall, the classification performance in GDMD spectrum was better than those in most of the remaining methods. In the statistical evaluation, there were significant differences in all classification performance and tasks ( > . × − p 3 63 10 8 ) using analysis of variance (ANOVA) or Friedman test. In the following evaluations, we indicate the post-hoc comparison results. Accuracy in GDMD spectrum was significantly higher than that of other methods ( > .
r 0 39 and < . p 0 03) for both recognition tasks, with the exception of GDMD Laplacian and DMD spectrum in the team-defence recognition ( > . r 0 16 and > .
p 0 05). AUC and F-measure in GDMD spectrum were also higher than those of the other methods ( > .
r 0 40 and < . p 0 028) for both recognition tasks, with the exception of hand-crafted feature in the team-defence recognition (both AUC and F-measure: > . r 0 34 and > . p 0 05). As a feature extraction method, Graph DMD was a better method than most of the other methods, and especially simple Graph DMD spectrum was better than the graph features in the collective motions.
Analysis of visualised Graph DMD modes. Since our method can decompose the data into temporal dynamics and spatial coherent structures in the observed data space, we can interpret the decomposed modes that contribute to the classification results. Here we show the averaged spectra for each label and classification task to identify the trends (note that the location of collective motion networks shown in Fig. 2 change with time; thus for an accurate understanding, we averaged and visualised them in an adjacency matrix form).
For the team-defence recognition task, a strong spectrum for zone defence in Fig. 4a was observed in the relationship among the nearest defenders (expanded in the figure) but it was not observed in person-to-person defence Fig. 4b. The variability of the feature vectors is also presented via boxplots in Fig. 5a-c. In examining the individual contribution, statistical analysis using logistic regression indicates that the relationships among defenders (D1-D2 and D4-D5: odds ratio > 31.81 and < .
p 0 023; D1-D3 was also significant but too small), those among attackers and defenders (A1-D4, A1-D5, A2-D5, A4-D1, and A5-D5: odds ratio > 0.02 and < .  Figure 5. Boxplots of feature vectors as elements of Graph DMD modes. Feature vectors as elements of Graph DMD modes for label 1 (red) and 2 (green) regarding the feature vectors used in the team-defence (a-c) and offence (d and e) recognition tasks are shown (label 1: zone defence and offence with screen-play; label 2: person-to-person defence and offence without screen-play). For the team-defence recognition task, the elements of GDMD modes (a) among defenders, (b) those among attackers and defenders and (c) those among defenders and the ring were selected. For the team-offence recognition task, those (d) among attackers and (e) those among attackers and defenders were selected. Notations are same as Fig. 4 www.nature.com/scientificreports www.nature.com/scientificreports/ For the team-offence recognition task, a strong spectrum for the offence with screen-play in Fig. 4c was observed in the relationship between the nearest attackers; however, it was not observed for the offence without screen-play Fig. 4d (the boxplot is shown in Fig. 5d,e). Results of logistic regression indicate that the relationship in D1-D2 (odds ratio = 1.06 × 10 4 and = . × − p 5 36 10 5 ) and those among attackers and defenders (A1-D5, A4-D3 and A5-D5: odds ratio > 23.59 and < . p 0 009; A5-D3 was also significant but too small) significantly explained the labels of team-offence. For more details, see Supplementary Table 2 and Supplementary Fig. 2. With these methods, we can interpret the types of interaction between agents (and environments) that contributes to focusing dynamic properties (i.e. frequency in this case) with a physical meaning. In Discussion section, we discuss further semantic interpretation.

Discussion
The objective of this paper was to understand the network dynamics for complex collective motions via physically-interpretable classification using data-driven spectral analysis of GDSs called Graph DMD. In this section, we discuss the comparison with the conventional approaches, semantically interpret the visualised Graph DMD modes, describe the limitation and the future perspectives for our methods. We then present our conclusions.
Using ballgame data as an example, we classified the strategic collective motions in the team-defence and offence recognition tasks more successfully than most of the existing methods. The results of the comparison with graph features (GDMD Laplacian and GCN) suggest that the collective motions in this case can be classified with the node information (nearest to the ball) rather than the graph features. The results of other DMD methods (DMD spectrum and KDMD spectral kernel) suggest that both the reflecting graph structure and the use of appropriate sliding windows were important for recognition. Furthermore, in comparison with KDMD spectral kernel 40 (i.e. decomposition in a feature space), the proposed method has an advantage in physically and semantically interpreting the DMD modes in the observed data space. Regarding an end-to-end neural network approach (i.e. spatio-temporal GCN), the model for motion capture data (i.e. joint location data for a few people) must be customised. However, each agent is not physically connected in this case; thus, the positional relationship (i.e. constraint condition) dynamically changes in contrast to the motion capture data. In this paper, we used a simple setting with minimal correction; however, further adaptation to complex collective motions is required (e.g. more appropriate customisation of the network and adequate learning such as using a larger amount of data). For hand-crafted features, we compute simple features to demonstrate the validity of our methods in this paper. The use of a more customised framework 44 which includes the detection of a few related people might improve specific classification performance. However, in this paper, we proposed a more generalised classification framework for global dynamics of collective motions, which successfully performed two different recognition tasks (team-defence and offence) by selecting only the elements of the Graph DMD modes as feature vectors. Furthermore, our method can obtain the dynamical properties of the collective motion network.
The proposed method can obtain physically-interpretable dynamical properties of network dynamics and can perform classification even for complex collective motions. For example, in Fig. 5, we found the label-specific stronger spectra in two recognition tasks, in which the relationships among the nearest attackers or defenders, those among attackers and defenders far from the ball and those between defenders and the ring, which were oscillated in low-frequency modes. Semantically, in zone defence, defenders adjusted their positions more according to a teammate (Fig. 5a), the attacker with the ball (Fig. 5b: especially in the distant case from the ball) and less according to the geometric reference point (Fig. 5c) than the defenders in person-to-person defence. For team-offence with screen-play, attackers adjusted their positions more according to a teammate (Fig. 5d) and the distant defender from the ball (Fig. 5e) than the attackers in team-offence without screen-play. Generally, our approach can be applied to the analysis of complex global dynamics in groups of living organisms or artificial agents, which currently eludes mathematical formulation. The human group in a sport used in this study can be considered as an example in which communication is likely to be measured in a physical space 38 . Therefore, by measuring the outputs of communication and assuming that there are underlying dynamics behind the obtained data, our approach can handle various means of communication between agents. In practice, in various sports teams or other human communities, supervisors (e.g. experimenters, coaches and teachers) spend considerable amounts of time analysing the collective motions in their domain. Application of a system, such as the one presented here, can lead to the creation of useful plans that are currently derived only from their implicit experience.
However, there are some limitations to this study. One is the validation of whether our approach can extract true dynamical properties if used as an equation-free method, which cannot confirm the true dynamical properties (e.g. frequencies with growth/decay rate), as well as in the previous works 31,32 . Originally, DMD demonstrates its strength for the dynamical systems which can be mathematically defined 41,48 or of which solutions are empirically known 29,30 . We instead validated our approach using classification performance, a qualitative evaluation with specific knowledge of the sport domain (e.g. low-frequency band) and the reconstruction error (see Materials and Methods). However, a general quantitative validation method for the unknown dynamics is further required. Another is to reflect the more local interaction dynamics such as local competitive and cooperative play by the attackers and defenders 22,[49][50][51] , which can provide more practical information in the sport domain. Although the purpose of Graph DMD is to extract the underlying global dynamics of GDSs and we can obtain the interpretable local spectra in the DMD modes, there are other approaches for extracting the more specific local dynamics, such as using heuristic selection of Graph DMD modes and additional hand-crafted features (e.g. specific geometrical information). Even when using only players' location data such as in this study, such methods reflecting local competition can be applied to more practical application such as score prediction 38,40 and prediction of a player to obtain the ball after shot 52 .
In conclusion, we applied a data-driven spectral analysis called Graph DMD, which obtains physically-interpretable dynamical properties of network dynamics for the classification even for complex Scientific RepoRtS | (2020) 10:3005 | https://doi.org/10.1038/s41598-020-58064-w www.nature.com/scientificreports www.nature.com/scientificreports/ collective motions. We classified the motions in different global behaviours and discovered that, in addition to the physical properties, the contextual node information was critical for classification. Furthermore, we discovered the label-specific stronger spectra in the relationship between the nearest agents, which provided physical and semantic interpretation. Potential applications to other collective motions are related to the motions which are difficult to discretize strategic motions such as football and handball, or other collective motions such as pedestrians 24 or animal group behaviour 38,53 . Our approach contributes to the understanding of complex networks involving collective motions from the perspective of nonlinear GDSs.

Materials and Methods
Positional data in a ballgame. The positional data of players and the ball (25 frames per second) were obtained from actual men's Asian international level practical games held in 2015 and preprocessed by STATS SportVU system (Northbrook, IL, USA). We obtained the consent to use it for research. We analysed 220 min of play (in 4 days) in which the two teams scored 746 points (386 vs. 360). For each day, players performed one and a half games (i.e. 60 min) except for 1 day (only one game). The positional data contained the XY position of each player and the XYZ coordinates of the ball on the court (All-court: 28 × 15 m; half-court: 14 × 15 m). After the following data segmentation, we obtained 319 attack-segments.
Data segmentation. Prior to data segmentation, we used Butterworth low-pass filter at 2 Hz and a custom-made automatic individual play-detection system to detect shots using positional data similar to that used in previous studies 22,38,44 . We analysed an attack-segment defined as the period that begins when all players enter the attacking side of the court and end before a shot (we analysed only the attack-segment finishing with a shot). Then, the authors, who have experience in playing and coaching basketball, manually labelled all attack-segments into a zone defence, a person-to-person defence as the team-defence recognition task and offence with and without screen-plays as the team-offence recognition task.
Person-to-person defence is a team-defence strategy in which defenders basically guard their predetermined attackers. In contrast, zone defence is another team-defence strategy in which defenders basically guard their predetermined areas. Although zone and person-to-person defence are exclusive by definition, they are actually mixed according to different situations. For example, against an attacker's penetration with the ball to the ring, the defenders in person-to-person defence also guard the area near the ring; in contrast, against an attacker's shot, the defenders in zone defence also guard the attacker with the ball. Despite this ambiguity, people with the experience of games can distinguish the two team-defence strategies by observing the defensive motion without the ball.
Screen-play is the basic and minimal strategic cooperative play in basketball 44 , in which an attacker stands on the course of defence player like a'screen' and prevents the defender from defending another attacker in a legal way. We labelled attack-segments into an offence with (at least including a screen-play) and without screen-play. Note that we focused on the global network dynamics of collective motions, thus we did not segment and label screen-plays themselves (i.e. in a spatio-temporally accurate sense such as in the previous study 44 ).
We randomly created a validation dataset (30 attack-segments) for the selection of parameters regarding Graph DMD, and a test dataset (289 attack-segments) for the subsequent analyses.
Creating adjacency matrix series. Next, we compute the adjacency matrix sequence in the attack-segment for inputting to the subsequent graph DMD. Here, using the Gaussian kernel in Eq. (5), we converted the distances between all individuals and the geometric reference position (i.e. the ring) into weights of adjacency matrix series, which indicate 1 if close to each other and 0 if far from each other. In this way, we can obtain more stable time series than those directly using the distances. The numerator in the exponential in Eq. (5) represents the distance between all individuals and the ring. Since the order of i cannot be uniquely determined in general, we must define the order based on reasonable grounds. For example, the order of playing positions (e.g. guard, forward and centre) or jersey numbers does not provide a fair comparison among attack-segments because the players themselves frequently changes throughout attack-segments. A simple way to determine the order is based on the distance between players and a reference point (e.g. ball or ring). However, if the order is simply determined by proximity to the reference point, there may arise a problem that, for example, the distance between the third and fifth players is shorter than that between those (the third and fifth) and the fourth players. Thus, first, we simply set = i 1, 2 for the two attackers in order from the attacker closer to the ball. Then, we set = i 3 for the closest attacker in terms of the sum of the distance to the above two attackers. In this way, we set = i 4, 5 for the closest attacker based on the sum of the distance to the three and four attackers, respectively. Similarly, we set = … i 6, , 10 for the five defenders and set = i 11 for the ring position to take the geometric information into consideration.
Next, denominator 2σ in the exponential in Eq. (5) is a coefficient for adjusting the value of A i j t , , . In this paper, we set σ for the relationship between players that satisfies = .
, 2 m (i.e. σ = . 1 5 2 /2 log 2), which is considered to be a meaningful proximity between players (e.g. 1 m is near and 2 m is far in terms of the contact of two players). For the additional analysis regarding σ and the ring information, see Supplementary Text 3. Also, we set σ for the relation between players and the ring that satisfies = .
11, 2 (i.e. σ = 6 2 /2 log 2). In this way, a sequence of the matrix A t with A i j t , , as an element at time t is created as a sequence of adjacent matrices of the graph. In this paper, we used this adjacency matrix series for both team-defence and offence recognition tasks.
Selection of an appropriate representation of the data is a fundamental problem in pattern recognition. Time-series data is challenging to design features for because of difficulties in reflecting the data structure (including time length). A simple method involves using the representative values (e.g. average and maximum values).

Basic DMD and its variants.
Here, we first describe the basic DMD procedure and briefly explain the variants of DMDs including Graph DMD. Among several possible methods to compute the spectral decomposition from data, DMD 29,30 is the most popular algorithm, which estimates an approximations of the decomposition in Eq. (2). Consider a finite-length observation sequence … τ y y y , , , . Then, DMD basically approximates it by calculating the eigendecomposition of matrix = † F YX , where X † is the pseudo-inverse of X. The matrix F may be intractable to analyse directly when the dimension is large. Therefore, in the popular implementation of DMD called exact DMD 47 , a rank-reduced representation F based on SVD is applied. That is, ( 1) , where * is the conjugate transpose. Thereafter, we perform eigendecomposition of  ∈ × F p p to obtain the set of the eigenvalues λˆj and eigenvectors w j . Then, we estimate the Koopman modes in Eq. (2): (see Results). Although A t is also vectorised, unlike SVD, M reflects order-3 tensor structure based on tensor-train decomposition. Therefore, both Koopman spectral analysis with dependent structure among observables and Graph DMD can be calculated without breaking the dependent structure.
Computation of Graph DMD. In this subsection, we describe the implementation of Graph DMD 33 . The procedure to compute Graph DMD modes is described in Results. The Matlab code we used is available at https:// github.com/keisuke198619/GraphDMD. To compute this, a modified tensor-based DMD 33,42 using tensor-train decomposition 43 is applied. The tensor-train decomposition is considered to be relatively stable and scalable for high-order tensors compared with the other tensor decomposition methods 43 . The basic algorithm of tensor-train decomposition for an order-N tensor (i.e.  ∈ ×…× A n n N 1 , where n l denotes the dimensionality of the l-th mode , by serial matricisations and SVDs. The input tensor of Graph DMD is an order-3 tensor; thus it is decomposed into a matrix, an order-3 tensor, and a matrix. For obtaining the time dynamics shown in Eq. (4) and Figs. 1 and 2, we need to compute b j,0 . Since here we do not need to consider the tensor structure in the time dynamics (i.e. they are univariate time-series), we compute b j,0 in a similar way of exact DMD by matricising the DMD mode = Z Z j j :,:, (such that the row is the time stamp) and vectorising the initial value A 0 . After computing the time dynamics, we can compute the reconstructed data shown in Eq. (4). Using the reconstructed data, we eliminated invalid sliding windows (described below) and investigated the validity of the decomposition. In this paper, as a criterion of valid decomposition, we used variability accounted for (VAF) which is commonly used in non-orthogonal dimensionality reduction 57 . VAF is defined as the square error of the reconstructed data and the original data such that where A and Â are matricised input data tensor A (such that the row is time) and its reconstructed tensor for the jthe mode, and ||·|| F is the Frobenius norm. In this paper, if the maximum VAF j < 0.01, the sliding window was eliminated from the subsequent analyses. We confirmed all attack-segments have at least one valid sliding window.
Next, we investigated the validity of the decomposition when changing the parameters (e.g. the Graph DMD tolerance, sliding window size and cutoff frequency) using a validation dataset. Here we describe the summary (for the details, see Supplementary Text 2). First, we set Graph DMD tolerance ε (i.e. the tolerance in the successive SVD in tensor-train decomposition) to 1.0 × 10 −5 based on the result in Supplementary Fig. 1a. For the sliding temporal windows, we set the window size to 50 frames (2 s) including overlaps of 25 frames (1 s) based on the result in Supplementary Fig. 1b. Regarding the cutoff frequency, we removed the DMD modes over 2 Hz based on the result in Supplementary Fig. 1c. Next, we describe the detailed computation of the feature vectors in classification using the Graph DMD modes.
Feature vectors using Graph DMD modes. We first adjusted the computed DMD modes. One is a normalisation within the DMD modes such that the maximum absolute value of the elements is 1 because the amplitude of the modes depends on that of the time dynamics. Second is symmetrisation. Since the computed Graph DMD modes were different from symmetric matrices in a precise sense (although the ideal modes are symmetric matrices), we added the modes to the transposed modes and divided them by two. Finally, to evaluate the global complex behaviours, we averaged the DMD modes for all valid sliding windows.
Classification. After Graph DMD, we perform classification in three different approaches. For details, see Supplementary Text 1. The first two approaches create feature vectors and the third automatically extracts the features via a neural network approach. The final two approaches computed graph features using existing methods: the second is graph Laplacian eigenvalues 39 as a baseline feature of a graph 45 (denoted GDMD Laplacian) and the third is deep graph convolutional neural networks 46 (denoted GDMD GCN) as a recent promising method for classifying the graph data. For all classification methods with the exception of the neural network approaches, we adopted logistic regression as a simple linear binary classification model.
As comparable methods for classifying graph sequence data, we adopted four methods. The first is a method using vectorised basic DMD modes 47 as a baseline of DMD approaches. The second is the Koopman spectral kernel 40 using DMD with reproducing kernels 41 as the existing method 38 for classifying the collective motion dynamics. The third is a hand-crafted feature as a simple baseline method. Fourth is an advanced neural network approach for classifying graph sequence data called spatio-temporal GCN 36 as a recent promising end-to-end method. For the detailed setups, see Supplementary Text 1.

Statistical analysis.
For an accurate evaluation of classification performance from the multiple viewpoints, in addition to accuracy, AUC and F-measure were calculated to compare the classification performance. AUC is based on the ROC curve, which is generated by plotting a cumulative distribution function of the true-positive rate with respect to the false-positive rate. Then, AUC takes various decision thresholds into consideration. Next, recall and precision rate were computed for F-measure. Recall rate is defined as the ratio of the sum of true positives and true negatives to the number of true positives (the true-positive rate), and the precision rate is defined as the ratio of the sum of true positives and true negatives to false positives. The trade-off curve between recall and precision was created using the cumulative distribution function. To evaluate the trade-off, the F-measure was calculated as follows: (2 × precision rate × recall rate)/(precision rate + recall rate).
For investigating the robustness of the classification results, we repeated test sessions for logistic regression five times using different test sets in analogy to five-fold cross-validation (i.e. classified 25 times). To compare the various methods, since the hypothesis of homogeneity of variances between methods was rejected with Levene's test, the Kruskal-Wallis nonparametric tests were performed. As the post-hoc comparison, Wilcoxon rank sum Scientific RepoRtS | (2020) 10:3005 | https://doi.org/10.1038/s41598-020-58064-w www.nature.com/scientificreports www.nature.com/scientificreports/ test with Bonferroni correction was used within the factor where a significant effect in Kruskal-Wallis test was found. Since comparisons were only performed between GDMD spectrum and others based on our hypothesis, p-value was multiplied by six. We used r values as the effect size for Wilcoxon rank sum test. We adopted logistic regression for simply investigating the effect of each element of the feature vector. We calculated the odds ratio (related to effect size), its p-value and 95% confidence interval. For all statistical calculations, < .
p 0 05 was considered significant. All statistical analyses were performed using the MATLAB 2018a Statistics and Machine Learning Toolbox (The MathWorks, Inc., MA, USA).