Identifying Oscillatory Hyperconnectivity and Hypoconnectivity Networks in Major Depression Using Coupled Tensor Decomposition

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Identifying Oscillatory Hyperconnectivity and Hypoconnectivity Networks in Major Depression Using Coupled Tensor Decomposition © 2021 the Authors Published version Liu, Wenya; Wang, Xiulin; Xu, Jing; Chang, Yi.; Hämäläinen, Timo; Cong, Fengyu Liu, W., Wang, X., Xu, J., Chang, Yi., Hämäläinen, T., & Cong, F. (2021). Identifying Oscillatory Hyperconnectivity and Hypoconnectivity Networks in Major Depression Using Coupled Tensor Decomposition. IEEE Transactions on Neural Systems and Rehabilitation Engineering, 29, 1895- 1904. https://doi.org/10.1109/tnsre.2021.3111564 2021 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 29, 2021 1895 Identifying Oscillatory Hyperconnectivity and Hypoconnectivity Networks in Major Depression Using Coupled Tensor Decomposition Wenya Liu , Xiulin Wang , Student Member, IEEE , Jing Xu, Yi Chang, Timo Hämäläinen , Senior Member, IEEE , and Fengyu Cong, Senior Member, IEEE Abstract —Previous researches demonstrate that major depression disorder (MDD) is associated with widespread network dysconnectivity, and the dynamics of functional connectivity networks are important to delineate the neural mechanisms of MDD. Neural oscillations exert a key role in coordinating the activity of remote brain regions, and various assemblies of oscillations can modulate different networks to support different cognitive tasks. Studies have demonstrated that the dysconnectivity of electroencephalography (EEG) oscillatory networks is related with MDD. In this study, we investigated the oscillatory hyperconnectivity and hypoconnectivity networks in MDD under Manuscript received January 11, 2021; revised May 19, 2021 and August 5, 2021; accepted September 7, 2021. Date of publication September 9, 2021; date of current version September 17, 2021. This work was supported in part by the National Natural Science Foundation of China under Grant 91748105, in part by the National Foundation in China under Grant JCKY2019110B009 and Grant 2020-JCJQ-JJ-252, in part by the Fundamental Research Funds for the Central Universities in Dalian University of Technology, China, under Grant DUT2019 and Grant DUT20LAB303, and in part by the Scholarships from China Scholarship Council under Grant 201706060263 and Grant 201706060262. (Corresponding author: Fengyu Cong.) This work involved human subjects or animals in its research. Approval of all ethical and experimental procedures and protocols was granted by the First Affiliated Hospital of Dalian Medical University. Wenya Liu is with the Faculty of Electronic Information and Electrical Engineering, School of Biomedical Engineering, Dalian University of Technology, Dalian 116024, China, and also with the Faculty of Information Technology, University of Jyväskylä, 40014 Jyväskylä, Finland (e-mail: [email protected]). Xiulin Wang is with the Department of Radiology, Affiliated Zhongshan Hospital of Dalian University, Dalian 116001, China, and also with the Faculty of Electronic Information and Electrical Engineering, School of Biomedical Engineering, Dalian University of Technology, Dalian 116024, China (e-mail: [email protected]). Jing Xu and Yi Chang are with the Department of Neurology and Psychiatry, First Affiliated Hospital, Dalian Medical University, Dalian 116011, China (e-mail: [email protected]; [email protected]). Timo Hämäläinen is with the Faculty of Information Technology, University of Jyväskylä, 40014 Jyväskylä, Finland (e-mail: [email protected]). Fengyu Cong is with the Faculty of Electronic Information and Electrical Engineering, School of Biomedical Engineering, Dalian University of Technology, Dalian 116024, China, also with the Faculty of Information Technology, University of Jyväskylä, 40014 Jyväskylä, Finland, also with the Faculty of Electronic Information and Electrical Engineering, School of Artificial Intelligence, Dalian University of Technology, Dalian 116024, China, and also with the Key Laboratory of Integrated Circuit and Biomedical Electronic System, Dalian University of Technology, Dalian, Liaoning 116024, China (e-mail: [email protected]). Digital Object Identifier 10.1109/TNSRE.2021.3111564 a naturalistic and continuous stimuli condition of music listening. With the assumption that the healthy group and the MDD group share similar brain topology from the same stimuli and also retain individual brain topology for group differences, we applied the coupled nonnegative tensor decomposition algorithm on two adjacency tensors with the dimension of time ×frequency ×connectivity ×subject, and imposed double-coupled constraints on spatial and spectral modes. The music-induced oscillatory networks were identified by a correlation analysis approach based on the permutation test between extracted temporal factors and musical features. We obtained three hyperconnectivity networks from the individual features of MDD and three hypoconnectivity networks from common features. The results demonstrated that the dysfunction of oscillatory networks could affect the involvement in music perception for MDD patients. Those oscillatory dysconnectivity networks may provide promising references to reveal the pathoconnectomicsof MDD and potential biomarkers for the diagnosis of MDD. Index Terms —Dynamic functional connectivity, coupled tensor decomposition, major depression disorder, naturalistic music stimuli, oscillatory networks. I. INTRODUCTION MAJOR depression disorder (MDD) is a globally common psychiatric disorder characterized by deficits of affective and cognitive functions [1]–[3]. It is almost a consensus to researchers that MDD is accompanied by abnormal functional connectivity (FC) between some brain regions, like cortical regions in the default mode network (DMN), rather than the aberrant response of individual brain regions [3]–[6]. Music therapy is associated with improvements in mood, which has made it an attractive tool for MDD treatment [7]. Previous studies have suggested that the oscillatory asymmetry and dysconnectivity could be the potential biomarkers of MDD during music perception [8]–[10]. An increasing amount of researches have demonstrated that FC presents the potential of temporal variability across different time-scales (from milliseconds to minutes) to support continuous cognitive tasks. This is termed as dynamic functional connectivity (dFC), and it represents the processes by which networks and subnetworks coalesce and dissolve over time, or cross-talk between networks [11]–[13]. Recently, researches have reported abnormal dFC of specific brain regions and This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 1896 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 29, 2021 neural networks in MDD using resting-state functional Magnetic Resonance Imaging (RS-fMRI) [3], [5], [13], [14]. For example, Demirtas et al. found a decreased variability of FC in the connections between the DMN and the frontoparietal network [5]. Kaiser et al. showed that MDD patients presented decreased dFC between medial prefrontal cortical (MPFC) regions and regions of parahippocampal gyrus within the DMN, but increased dFC between MPFC regions and regions of insula. They showed that MDD was related to abnormal patterns of fluctuating communication among brain systems involved in regulating attention and self-referential thinking [13]. The decreased dFC variability was reported between anterior DMN and right central executive network (CEN) in MDD, which indicated a decreased information processing and communication ability [14]. Existing researches about dFC in MDD mostly focus on resting-state conditions. However, little is known about the abnormalities of dFC during music listening conditions. Benefiting from the high temporal resolution, electroencephalography (EEG) can record electrical brain activity dynamics at a millisecond scale with rich frequency contents. The oscillation acts as a bridge to connect different brain regions with resonant communication, which can regulate changes of neuronal networks and cause qualitative transitions between modes of information processing [15]–[17]. Impaired coordination of brain activity associated with abnormal electrophysiological oscillations contributes to the generation of psychaitric disorders [18]. Numerous studies have investigated EEG oscillatory FC of MDD in resting-state, and dysconnectivity networks are mostly notable in theta, alpha and beta oscillations [16], [19], [20]. However, most previous studies filter EEG signals into a range of frequency bands (e.g., 8-13 Hz for the alpha band), and ignore the exhaustive spectral dynamics in FC [19], [20]. Music perception is a complex cognitive task, which is characterized with dynamics of frequency-specific brain networks for musical features processing [21]–[25]. To the best of our knowledge, the oscillatory dFC in MDD during music perception has not been well investigated yet. Considering the temporal dynamics and spectral modulations of spatial couplings (e.g., functional connectivity) for multiple participants in a cognitive task, a multi-way dataset structure is naturally formed. This multi-dimensional nature points to the adoption of tensor decomposition models instead of matrix decomposition models, which normally fold some dimensions and ignore the hidden interactions across different modes [24], [26]–[29]. Canonical Polyadic (CP) decomposition is derived in terms of the sum of multiple rank-one tensors, and each rank-one tensor represents the covariation of the corresponding components from each mode [30], [31]. The CP decomposition is well implemented into the extraction of multi-mode EEG features from the multiway dataset (e.g., channel ×frequency ×time ×subject) [31]–[34]. Recently, Zhu et al. applied CP decomposition to explore the task-related dFC characterized by spatio-temporal- spectral modes of covariation from the adjacency tensor (connectivity ×time-subject ×frequency) [23], [35]. However, those applications only focus on the decomposition of one single tensor, which are based on the assumption that the underlying spatio-spectral features are consistent among subjects or groups [25], [29]. Coupled tensor decomposition (CTD), the extension of tensor decomposition to multiple block tensors, enables the simultaneous extraction of common features shared among tensors and individual features specified for each tensor. For biomedical data, the coupled matrix, matrix-tensor or tensor decomposition (also known as linked component analysis) are mostly used for data fusion [36]–[38]. However, to the best of our knowledge, no studies have used CTD to investigate the pathologic networks of MDD or other psychiatric disorders. In our study, we applied a low-rank double-coupled nonnegative tensor decomposition (DC-NTD) model to explore the temporal and spectral dynamics of spatial couplings in MDD during music listening. The proposed analysis pipeline is totally data-driven. We analyzed the whole-brain FC to avoid prior knowledge about regions of interest, and we investigated the exhaustive assemblies of oscillations to avoid the selection of the frequency range. Figure 1 shows the diagram of the analysis pipeline of this study. In this paper, scalars, vectors, matrices and tensors are denoted by lowercase, boldface lowercase, boldface uppercase and boldface script letters, respectively, e.g., x,x, X,X. Indices range from 1 to their capital version, e.g., i=1,··· ,I. II. MATERIALS AND METHODS A. Simulated Data To validate the feasibility of the proposed method, we firstly applied it on the simulated data. Two tensors with the size of 500 ×59 ×2278, representing time ×frequency × connectivity, were created as follows: X1=˜ X1+N1 =t1◦f1◦c1+t2◦f2◦c2+t3◦f3◦c3+N1,(1) X2=˜ X2+N2 =t4◦f1◦c1+t5◦f2◦c2+t6◦f4◦c4+N2,(2) where ˜ Xm,m=1,2 represented the ground truth networks, and Nn,n=1,2 were the nonnegative noise created by the absolute values of white noise with the size of 500 × 59 ×2278. In the time domain, each temporal component ti,i=1,2,··· ,6 was simulated by the absolute value of white noise to ensure the nonnegativity of the synthetic tensor X, and no coupled temporal component existed between two tensors. In the frequency domain, four spectral components fj,j=1,2,··· ,4 were constructed by Hanning windows and white noise with bandwidth centered at 5 Hz, 10 Hz, 15 Hz and 20 Hz, and two spectral components were set to be coupled between two tensors. In the adjacency domain, four adjacency components ck,k=1,2,··· ,4, representing auditory network (AUD), visual network (VIS), salience network (SAN), and dorsal attentional network (DAN), were constructed with the Desikan-Killiany anatomical atlas according to Kabbara’s work [39], and two adjacency components were coupled between two tensors. The synthetic data were shown in Figure 2(a). LIU et al. : IDENTIFYING OSCILLATORY HYPERCONNECTIVITY AND HYPOCONNECTIVITY NETWORKS 1897 Fig. 1. Diagram of the analysis pipeline. (a) Adjcency matrix construction in each time window and each frequency bin. After source reconstruction, the cortical signals were segmented by overlapping time windows, and wavelet transform was applied for each time course within each time window. Phase lag index was used to obtain the adjacency matrix for each time window and each frequency bin. (b) Adjacency tensor construction and decomposition. A 4-D adjacency tensor was constructed for each group with the dimension of time ×frequency ×connectivity ×subject, and coupled tensor decomposition was implemented with coupled constraints in spectral and adjacency modes. The 4-D core tensor is superdiagonal with values of 1. (c) The identification of hyperconnectivity and hypoconnectivity networks by music modulation. Five musical features were extracted with MIR toolbox from tango music, and correlation analysis was conducted between musical features and decomposed temporal factors to identify music-induced brain networks. Hyperconnectivity and hypoconnectivity networks were summarized from the results of music modulation. B. EEG Data Description 1) Participants: Twenty MDD patients and nineteen healthy controls (HC) participated in this experiment. All the patients were from the First Affiliated Hospital of Dalian Medical University in China. This study has been approved by the ethics committee of the hospital, and all participants signed the informed consent before their enrollment. None of the participants has reported hearing loss and formal training in music. All the MDD patients were primarily diagnosed by a clinical expert and tested according to Hamilton Rating Scale for Depression (HRSD), Hamilton Anxiety Rating Scale (HAMA) and Mini-Mental State Examination (MMSE). The means and standard deviations (SD) of age, education, clinical measures, duration of illness and gender for both groups were listed in Table I. 2) EEG Data: During the experiment, participants were told to sit comfortably in the chair and listen to a piece of music. A 512-second musical piece of modern tango Adios Nonino by Astor Piazzolla was used as the stimulus due to its rich musical structure and high range of variation in musical features such as dynamics, timbre, tonality and rhythm [21], [40]. The EEG data were recorded by the Neuroscan Quik-cap device with 64 electrodes arranged according to the international 10-20 system at the sampling frequency of 1000 Hz. The electrodes placed at the left and right earlobes were used as the references. TABLE I MEANS AND STANDARD DEVIATIONS OF AGE,EDUCATION,CLINICAL MEASURES,DURATION OF ILLNESS AND GENDER FOR THE HC GROUP AND THE MDD GROUP The data were visually checked to remove obvious artifacts from head movement and down-sampled to fs =256 Hz for further processing. Then 50 Hz notch filter and highpass and low-pass filters with 1 Hz and 30 Hz cutoff were applied. We interpolated the bad intervals of one channel by 1898 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 29, 2021 the mean values of their spherical adjacent channels. Eye movements artifacts were rejected by independent component analysis (ICA). 3) Musical Features: In this study, two tonal and three rhythmic features were extracted by a frame-by-frame analysis approach using MIR toolbox [41]. The duration of each frame was 3 seconds, and the overlap between two adjacency frames was 2 seconds. Therefore, we got 510 samples for the time courses of each musical feature at a sampling frequency of 1 Hz. In this study, we only used the first T=500 samples of each musical feature due to the length of recorded EEG data. Tonal features include Mode and Key Clarity, which represent the strength of major of minor mode and the measure of tonal clarity, respectively. Rhythmic features include Fluctuation Centroid, Fluctuation Entropy, and Pulse Clarity. Fluctuation centroid is the geometric mean of the fluctuation spectrum representing the global repartition of rhythm periodicities within the range of 0–10 Hz. Fluctuation Entropy is the Shannon entropy of the fluctuation spectrum representing the global repartition of rhythm periodicities. Pulse Clarity is an estimate of clarity of the pulse. C. Source Reconstruction Source reconstruction procedure was performed with opensource Brainstorm software [42]. For forward modeling, we used the symmetric boundary element method (BEM) to compute the volume-conductor model with the MNI-ICBM152 template corresponding to a grid of 15000 cortical sources. For source modeling, minimum norm estimate (MNE) was applied with a measure of the current density map and constrained dipole orientations (normal to cortex). Then, the Desikan- Killiany anatomical atlas was used to parcellate the cortical surface into C=68 regions, and the principal component analysis (PCA) method was performed to construct the time course for each brain region. D. Dynamic Functional Connectivity Many studies have reported that the communication of brain regions or neural populations depends on phase interactions for electrophysiological neuroimaging techniques, like EEG [43]. To avoid source leakage, the pairwise synchronization was estimated by PLI to map the whole-brain FC [44]. In this study, to assess the dFC across both time and frequency, we segmented the source-space data into W=500 windows by the sliding window technique with a window length of 3 s and an overlap of 2 s according to the extraction framework of musical features. Then, we computed the time-frequency (TF) decomposition within each time window by the continuous wavelet transform with Morlet wavelets as basis function . We set the frequency bins as 0.5 Hz, and obtained F=59 samples in frequency domain in the range of 1-30 Hz. For the time window w, we can get the complex TF representation Pw∈RTw×Ffrom wavelet transform, where Tw=3fs, and the time and frequency-dependent phase at time twand frequency fcan be obtained by ϕ(tw,f)=arctan imag(Pw(tw,f)) real(Pw(tw,f)) ,(3) where imag() and real() represent the imaginary part and the real part of a complex value, respectively. For brain regions i and j, PLI can be computed as PLIi,j(w, f)=1 Tw       Tw  tw=1 sign(ϕi,j(tw,f))      ,(4) where ϕi,j(tw,f)) =ϕi(tw,f)−ϕj(tw,f)is the phase difference of brain regions iand jat time twand frequency f in time window w. Therefore, for each time window and each frequency point, we can form an adjacency matrix A∈RC×C, where Cmeans the number of brain regions. Because of the symmetry of FC matrix, we took the upper triangle of Aand vectorized it to a∈RN×1,whereN=C(C−1)/2=2278 represents the number of unique connections. Then, we can construct two adjacency tensors with the dimension of time × frequency ×connectivity ×subject, XHC ∈RW×F×N×SHC (500 ×59 ×2278 ×19) for the HC group and XMDD ∈ RW×F×N×SMDD (500 ×59 ×2278 ×20) for the MDD group, where SHC =19 and SMDD =20 mean the number of subjects in the HC group and the MDD group, respectively. E. The Application of Low-Rank Coupled Tensor Decomposition Considering the high computation load, the nonnegativity of the tensors (constrained to [0,1] due to PLI index) and high correlations in spatial and spectral modes, we applied a lowrank DC-NTD model which was more flexible to add desired constraints. 1) Low-Rank Coupled Tensor Decomposition: With the constructed tensors XHC ∈RW×F×N×SHC and XMDD ∈ RW×F×N×SMDD , the corresponding CP decomposition can be represented as XHC ≃RHC r=1u(1) r◦u(2) r◦u(3) r◦u(4) r= [[ U(1),U(2),U(3),U(4)]] and XMDD ≃RMDD r=1v(1) r◦v(2) r◦ v(3) r◦v(4) r=[[V(1),V(2),V(3),V(4)]] ,where◦denotes the vector outer product. u(n) rand v(n) rdenote the rth component of factor matrices U(n)and V(n),n=1,2,3,4, in the modes of time, frequency, connectivity and subject for two groups. RHC and RMDD are the ranks of XHC and XMDD, respectively. Considering the nonnegativity of constructed tensors and the coupled constraints in spectral and adjacency modes, we formulate it as a double-coupled nonnegative tensor decomposition (DC-NTD) model, where XHC and XMDD can be jointly analyzed by minimizing the following objective function: J(u(n) r,v(n) r) =XHC − RHC  r=1 u(1) r◦u(2) r◦u(3) r◦u(4) r2 F +XMDD − RMDD  r=1 v(1) r◦v(2) r◦v(3) r◦v(4) r2 F s.t.u(2) r=v(2) r(r≤Lf), u(3) r=v(3) r(r≤Lc). (5) · Fdenotes the Frobenius norm. Lfand Lcdenote the number of components coupled in spectral and adjacency modes, and Lf,c≤min(RHC,RMDD). The fast hierarchical LIU et al. : IDENTIFYING OSCILLATORY HYPERCONNECTIVITY AND HYPOCONNECTIVITY NETWORKS 1899 alternative least squares (FHALS), an accelerated version of the hierarchical alternative least squares (HALS) algorithm, has been effectively applied to a number of (coupled) tensor decomposition problems [25], [45], [46]. In this study, we apply the FHALS algorithm to optimize the DC-NTD problem in (5), and introduce the low-rank approximation to reduce computational complexity [47], [48]. Through the FHALS algorithm, the minimization problem in (5) can be converted into max(RHC,RMDD)rank-1 subproblems, which can be solved sequentially and iteratively. We take the rth subproblem as an example: min Jr=YHC r−u(1) r◦u(2) r◦u(3) r◦u(4) r2 F +YMDD r−v(1) r◦v(2) r◦v(3) r◦v(4) r2 F,(6) where YHC r . =XHC−RHC k=ru(1) k◦u(2) k◦u(3) k◦u(4) kand YMDD r . = XMDD −RMDD k=rv(1) k◦v(2) k◦v(3) k◦v(4) k. When calculating one of the variables, we need to fix the other variables and let the corresponding derivative be zero. For example, to determine u(n) r,welet∂Jr/∂u(n) rbe zero, and then we can obtain the following solution: u(n) r=YHC r,(n)ur−n/urTur−n,(7) where YHC r,(n)is the mode-nmatricization of YHC r.and  denote the Khatri-Rao product and Hadamard (element-wise) product. ur−n=u(4) r···u(n+1 r)u(n−1 r)···u(1) rand urTur−n=[ur]−nT[ur]−n.TakingYHC r,(n)=XHC (n)− U(n)[U−n]T+u(n) r[u−n r]Tinto (7), we can get u(n) r=u(n) r+XHC (n)u−n r−U(n)(n) r/(n) (r,r),(8) where XHC (n)is the mode-nmatricization of XHC and (n)=UTU−n. Suppose that the rank- ˜ RHC approximation of XHC obtained by unconstrained tensor factorization is expressed as [[ ˜ U(1),˜ U(2),˜ U(3),˜ U(4)]] ,˜ RHC ≤RHC, thus the mode-nunfolding of XHC can be written as XHC (n)= ˜ U(n)˜ U−nT . Therefore, the learning rule of u(n) rcan be reformulated as follows: u(n) r=u(n) r+˜ U(n)˜ (n) r−U(n)(n) r/(r,r),(9) where ˜ (n)=[ ˜ UTU]−n. Analogously, we can obtain the learning rule of v(n) ras follows: v(n) r=v(n) r+˜ V(n)˜ (n) r−V(n)(n) r/(n) (r,r),(10) where (n)=VTV−nand ˜ (n)=[ ˜ VTV]−n. The rank- ˜ RMDD approximation of XMDD is expressed as [[ ˜ V(1),˜ V(2),˜ V(3),˜ V(4)]] ,˜ RMDD ≤RMDD. Specially, if r≤ Lf,u(2) r=v(2) rand if r≤Lc,u(3) r=v(3) r, thus their solutions should be calculated as: u(2) r=v(2) r =u(2) r+˜ U(2)˜ (2) r−U(2)(2) r+˜ V(2)˜ (2) r−V(2)(2) r /(2) (r,r)+(2) (r,r),(11) and u(3) r=v(3) r =u(3) r+˜ U(3)˜ (3) r−U(3)(3) r+˜ V(3)˜ (3) r−V(3)(3) r /(3) (r,r)+(3) (r,r),(12) In order to obtain the nonnegative components, a simple “half-rectifying” nonlinear projection is applied. We update u(n) rand v(n) rsuccessively in each subproblem, and the max(RHC,RMDD)subproblems are optimized alternatively one after another until convergence. In this study, we adopt alternating least squares (ALS, [49]) algorithm to perform lowrank approximation. The FHALS-based DC-NTD algorithm is summarized in Algorithm 1. Algorithm 1 DC-NTD-FHALS Algorithm Input:XHC,XMDD,Lf,Lc,RHC,RMDD,˜ RHC,˜ RMDD 1Initialization: U(n),V(n),n=1,2,3,4 2Calculate ˜ U(n),˜ V(n),n=1,2,3,4 via unconstrained ALS 3while unconvergence do 4for n=1,2,··· ,4do 5for r=1,2,···max(RHC,RMDD)do 6Update u(n) r,u(n) rvia (9), (10), (11) and (12) 7end 8end 9end Output:U(n),V(n),n=1,2,3,4 2) Selection of Components: In this section, we will describe how to determine the number of totally extracted components RHC and RMDD, which refers to the hidden information in low-dimensional space for each block data, and the number of coupled components Lfand Lc, which reveal the common features between two-block data. For the selection of RHC and RMDD, we performed PCA on the matricization data X(3)∈ RF×WNS unfolded along frequency mode for each block data, and kept the number of components with 95% explained variance. The selection of the number of coupled components is a key issue for the conduction of the DC-NTD-FHALS algorithm and the explanation of results, and it always becomes an open issue depending on practical applications. In this study, we performed the fourth-order CP tensor decomposition based on the FHALS algorithm on two-block data separately, and calculated the correlation maps of extracted components between two-block data in spectral and adjacency modes, respectively. According to the correlation maps, we will select the number of highly correlated (coupled) components. The detailed implication procedure will be described in the results section. F. Identification of Music-Induced Hyperconnectivity and Hypoconnectivity Networks After the conduction of the low-rank DC-NTD-FHALS algorithm, we need to identify the music-induced oscillatory 1900 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 29, 2021 Fig. 2. Simulation illustration. (a) Three spatio-temporal-spectral patterns were simulated for two groups, and the first two patterns were identical in adjacency and spectral modes. (b) The reconstructed spatio-temporal-spectral patterns. networks that are abnormal involved in the MDD group. We conducted a correlation analysis approach between temporal factors and five musical features with Pearson correlation based on the permutation test method. To ensure the statistical significance of the correlation and consider the problem of multiple comparison, the Monte Carlo method and permutation test were applied to compute a significant threshold of correlation for each musical feature [21]–[23], [25]. For the time course of each musical feature, we kept the real part and replaced the imaginary part with random uniformly distributed phases, and performed Pearson correlation with the time courses of the extracted temporal components. Then, we repeated this procedure 100000 times, and obtained the threshold for each musical feature at a significant level of pcorrected <0.05. The coupled spectral and adjacency components are common features between CON and MDD groups, and the remaining components are individual features of each group. The oscillatory networks among common features that are involved in music perception in the HC group but not in the MDD group are identified as hypoconnectivity networks, and the oscillatory networks among individual features that are involved in music perception in the MDD group are identified as hyperconnectivity. III. RESULTS A. Results of Simulated Data We implemented the low-rank DC-NTD-FHALS algorithm on the synthetic data, and we set SNR =15, Lf=Lc=2, Fig. 3. Correlation analysis of the spectral factor and the adjacency factor extracted from the 4-D tensor decomposition for each block data. and RHC =RMDD =3. The extracted temporal, spectral and adjacency factors were shown in Figure 2(b). We ran 10 times of the low-rank DC-NTD-FHALS algorithm, and we obtained stable decomposition results with an averaged tensor fit of 0.864 and an averaged running time of 113.27 seconds. B. Results of EEG Data Through PCA analysis on the unfolded data along the spectral mode for two-block data, we extracted RHC =RMDD = 27 components for both HC group and MDD group. Then we performed the fourth-order CP tensor decomposition on each block data, and computed the correlation maps of spectral and adjacency factors between two groups, as shown in Figure 3. According to the correlation maps, we set the number of coupled components in the spectral mode Lf=25 and the number of coupled components in the adjacency mode Lc=7. We ran 20 times of the low-rank DC-NTF-FHALS algorithm, LIU et al. : IDENTIFYING OSCILLATORY HYPERCONNECTIVITY AND HYPOCONNECTIVITY NETWORKS 1901 Fig. 4. Three oscillatory hyperconnectivity networks. (a) Adjacency matrix representation of the network. The 68 brain regions are ordered from the left hemisphere to the right hemisphere. Within each hemisphere, the brain regions are arranged in the order of frontal lobe, temporal lobe, parietal lobe, and occipital lobe, as indicated in red, yellow, green, and blue color, respectively. Within each lobe, the brain regions are ordered according to their y-location from anterior regions to posterior regions. (b) The spectral component of the network. (c) Cortical space representation of the network in Lateral, medial and dorsal view. The networks I,II, and III are related to the musical features of Fluctuation Centroid, Fluctuation Centroid, and Key Clarity, respectively. and the averaged running time was 12616 seconds. The running time was 63819 seconds by one implementation of the DC-NTF-FHALS algorithm without the low-rank approximation, which indicated that the low-rank approximation could greatly reduce the computational load. After applying the low-rank DC-NTF-FHALS algorithm and correlation analysis with musical features, we summarized the results of 20 times of implementation and obtained three oscillatory hyperconnectivity networks, as shown in Figure 4, and three oscillatory hypoconnectivity networks, as shown in Figure 5. For hyperconnectivity networks, Figure 4I shows a right hemisphere dominated network modulated by oscillations of alpha and beta (10-16 Hz) bands and the musical feature of Fluctuation Centroid. The strong connections of this network connect the core regions of DMN, including medial prefrontal cortex (mPFC), precuneus cortex, and posterior cingulate cortex (PCC). Figure 4II indicates a left auditoryrelated network modulated by delta oscillations and the Fluctuation Centroid feature. An aberrant delta-specific prefrontal network is identified, which is related to the musical feature of Key Clarity, as shown in Figure 4III. For hypoconnectivity networks, Figure 5I and Figure 5II exhibit fronto-parietal networks which are mainly related to attention control. The fronto-parietal networks are modulated by oscillations of 8- 14 Hz and 10-19 Hz and musical features of Mode and Fluctuation Entropy, respectively. Figure 5III shows a lowfrequency (delta oscillations) modulated prefrontal network which is significantly related to the musical feature of Mode, and this network is implicated in complex cognitive functions. IV. DISCUSSION As far as we know, this study is the first attempt to investigate the aberrant dFC across temporal evolution and spectral modulation in MDD during music listening based on a coupled tensor decomposition approach. This study proposed a comprehensive framework to extract the FC networks characterized by spatio-temporal-spectral modes of covariation. We summarized three overactive oscillatory networks and three underactive oscillatory networks according to the analysis of musical modulations. MDD is characterized with imbalanced communications among large-scale functional networks, including hyperconnectivity and hypoconnectivity within specific brain networks or between distinct brain networks during restingstate, see a meta-analysis in study [6]. In our study, we also found hyperconnectivity and hypoconnectivity functional networks in naturalistic music perception. We identified a right hemisphere dominated hyperconnectivity network which involved the essential regions of DMN, including mPFC, PCC and precuneus cortex, as shown in Figure 4I. 1902 IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING, VOL. 29, 2021 Fig. 5. Three oscillatory hypoconnectivity networks. (a) Adjacency matrix representation of the network. The 68 brain regions are ordered from left hemisphere to right hemisphere. Within each hemisphere, the brain regions are arranged in the order of frontal lobe, temporal lobe, parietal lobe, and occipital lobe, as indicated in red, yellow, green, and blue color, respectively. Within each lobe, the brain regions are ordered according to their y-location from anterior regions to posterior regions. (b) The spectral component of the network. (c) Cortical space representation of the network in Lateral, medial and dorsal view. The networks I,II, and III are related to the musical features of Mode, Fluctuation Entropy, and Fluctuation Entropy, respectively. The hyperconnectivity in DMN are often considered as reflecting rumination, where MDD patients perseverate on negative, self-referential thoughts [50], [51]. Many researches have reported the hyperconnectivity within DMN in MDD, which supports that within-DMN hyperconnectivity is related to enhanced the positive connectivity in MDD [6], [51]. Figure 4II shows a delta band-modulated and left auditoryrelated network, which is activated by a rhythmic feature of Fluctuation Centroid. The delta band was demonstrated to have a substantial influence on the identification of natural speech fragments in a MEG study [52], and the decoding of rhythmic features was found to be significantly correlated with the auditory cortex during music perception [21], [53]. The abnormal delta band-modulated and left auditory-related network identified in our study might indicate that MDD patients were less involved in music perception. We identified two delta band modulated prefrontal networks, both of which were related to tonal features, Key Clarity and Mode, as shown in Figure 4III and Figure 5III. However, the prefrontal network in Figure 4III was hyperactive and right hemisphere lateralized, and the prefrontal network in Figure 5III was hypoactive and left hemisphere lateralized. The prefrontal cortex has been implicated in planning complex cognitive behavior, decision making and working memory. There are numerous lines of evidence demonstrating that prefrontal cortex is dysregulated in depression, and both increased and decreased functional connections in the prefrontal network may lead to the failure of inhibitory control in depression [54]–[57]. Those two prefrontal networks also have abnormal connections with temporal poles, which may indicate the dysfunction in semantic integration during music listening [58], [59]. Our findings are well supported by those literatures that the dysconnectivity in the prefrontal network can influence the high-order cognitive functions and information integration during music perception in MDD. Figure 5I and Figure 5II indicate hypoconnectivity fronto-parietal networks modulated by different oscillations and musical features. The abnormal development of the fronto-parietal network is a common feature across many psychiatric disorders with the deficit in cognitive control. Previous studies have demonstrated that MDD is characterized by hypoconnectivity within the frontoparietal network, which is involved in the top-down modulation of attention and emotion [6], [19], [60]. In our study, the key issue of applying coupled tensor decomposition is the selection of the number of all the extracted components and the number of coupled components. There are several methods for the selection of the number of extracted components in tensor/matrix decomposition, such as