This is a self-archived version of an original arcle. This version
may dier from the original in paginaon 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 differ-
ent networks to support different cognitive tasks. Studies
have demonstrated that the dysconnectivity of electroen-
cephalography (EEG) oscillatory networks is related with
MDD. In this study, we investigated the oscillatory hyper-
connectivity 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.
(Corre-
sponding 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 Infor-
mation Technology, University of Jyväskylä, 40014 Jyväskylä, Finland
(e-mail: wenyaliu0912@foxmail.com).
Xiulin Wang is with the Department of Radiology, Affiliated Zhong-
shan 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: xiulin.wang@foxmail.com).
Jing Xu and Yi Chang are with the Department of Neurology and
Psychiatry, First Affiliated Hospital, Dalian Medical Univer-
sity, Dalian 116011, China (e-mail: xujing.doc@aliyun.com;
changee99@gmail.com).
Timo Hämäläinen is with the Faculty of Information Technol-
ogy, University of Jyväskylä, 40014 Jyväskylä, Finland (e-mail:
timo.t.hamalainen@jyu.fi).
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: cong@dlut.edu.cn).
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 percep-
tion 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, natural-
istic music stimuli, oscillatory networks.
I. INTRODUCTION
M
AJOR depression disorder (MDD) is a g lobally com-
mon psychiatric disor der characterized by deficits of
affective and cognitive functions [1]–[3]. It is almost a con-
sensus 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 differ-
ent time-scales (from milliseconds to minutes) to support con-
tinuous 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, r esearches
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 Mag-
netic 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 r esearches about dFC
in MDD mostly focus on resting-state conditions. However,
little is known about the abnormalities o f dFC during music
listening conditions.
Benefiting from the high temporal resolution, electroen-
cephalography (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 neuron al networks and cause qua litative transitions
between modes of information processing [15]–[17]. Impaired
coordination of brain activity associated with abnormal elec-
trophysiological oscillations contributes to the generation of
psychaitric disorders [18]. Numerous studies have investigated
EEG oscillatory FC of MDD in r esting-state, and d yscon-
nectivity 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 exhaus-
tive 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 modula-
tions 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 decom position 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) decompo-
sition 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 extrac-
tion 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]. How-
ever, those applications only focus on the decomp osition 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 b est 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 non-
negative 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 d a ta-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 upper-
case 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. M
ATERIALS 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:
X
1
=
˜
X
1
+ N
1
= t
1
◦ f
1
◦ c
1
+ t
2
◦ f
2
◦ c
2
+ t
3
◦ f
3
◦ c
3
+ N
1
, (1)
X
2
=
˜
X
2
+ N
2
= t
4
◦ f
1
◦ c
1
+ t
5
◦ f
2
◦ c
2
+ t
6
◦ f
4
◦ c
4
+ N
2
, (2)
where
˜
X
m
, m = 1, 2 represented the ground truth networks,
and N
n
, 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
t
i
, 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
f
j
, j = 1, 2, ··· , 4 were constructed by Hanning windows
and white noise with ban dwidth centere d 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 c
k
, k = 1, 2, ··· , 4, representing
auditory network (AUD), visual network (VIS), salience net-
work (SAN), and dorsal attentional network (DAN), were con-
structed with th e 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 h ospital, 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 o f 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
M
EANS AND STA NDAR D DEVIATIONS OF AGE,EDUCATION,CLINICAL
MEASURES,DURATION OF ILLNESS AND GENDER FOR THE
HC GROUP AND THE MDD GROUP
The data were v isually checked to remove obvious artifacts
from head movement and down-sampled to fs = 256 Hz
for further processing. Then 50 Hz notch filter and high-
pass and low-pass filters with 1 Hz and 30 Hz cuto ff were
applied. We in terpolated 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 o f 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 mea-
sure 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 period-
icities 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 open-
source 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 p airwise 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 P
w
∈ R
T
w
×F
from wavelet transform, where
T
w
= 3 fs, and the time and frequency-dependent phase at
time t
w
and frequency f can be obtained by
ϕ(t
w
, f ) = arctan
imag(P
w
(t
w
, f ))
real(P
w
(t
w
, 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
PLI
i, j
(w, f ) =
1
T
w
T
w
t
w
=1
si gn(ϕ
i, j
(t
w
, f ))
, (4)
where ϕ
i, j
(t
w
, f )) = ϕ
i
(t
w
, f ) − ϕ
j
(t
w
, f ) is the phase
difference of brain regions i and j at time t
w
and frequency f
in time window w. Therefore, for each time window and each
frequency point, we can form an adjacency matrix A ∈ R
C×C
,
where C means the number of brain regions. Because of the
symmetry of FC matrix, we took the upper triangle of A and
vectorized it to a ∈ R
N×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, X
HC
∈ R
W ×F ×N ×S
HC
(500 × 59 × 2278 × 19) for the HC group and X
MDD
∈
R
W ×F ×N ×S
MDD
(500 × 59 × 2278 × 20) for the MDD group,
where S
HC
= 19 and S
MDD
= 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 low-
rank DC-NTD model which was more flexible to add desired
constraints.
1) Low-Rank Coupled Tensor Decomposition:
With the con-
structed tensors X
HC
∈ R
W ×F ×N ×S
HC
and X
MDD
∈
R
W ×F ×N ×S
MDD
, the corresponding CP decomposition can be
represented as X
HC
R
HC
r=1
u
(1)
r
◦ u
(2)
r
◦ u
(3)
r
◦ u
(4)
r
=
[[ U
(1)
, U
(2)
, U
(3)
, U
(4)
]] and X
MDD
R
MDD
r=1
v
(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)
r
and v
(n)
r
denote the r th 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. R
HC
and R
MDD
are the ranks of X
HC
and X
MDD
,
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) mod e l, where X
HC
and X
MDD
can
be jointly analyzed by minimizing the following objective
function:
J (u
(n)
r
, v
(n)
r
)
=X
HC
−
R
HC
r=1
u
(1)
r
◦ u
(2)
r
◦ u
(3)
r
◦ u
(4)
r
2
F
+X
MDD
−
R
MDD
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 ≤ L
f
), u
(3)
r
= v
(3)
r
(r ≤ L
c
). (5)
·
F
denotes the Frobenius norm. L
f
and L
c
denote the
number of components coupled in spectral and adjacency
modes, and L
f,c
≤ min(R
HC
, R
MDD
). The fast hierarchical
