Further author information: (Send correspondence to Ali Vosoughi)
Ali Vosoughi: E-mail<EMAIL_ADDRESS>
# Large-scale Augmented Granger Causality (lsAGC) for Connectivity Analysis in
Complex Systems: From Computer Simulations to Functional MRI (fMRI)

Axel Wismüller  Department of Electrical and Computer Engineering, University
of Rochester, NY, USA  Department of Imaging Sciences, University of
Rochester, NY, USA  Department of Biomedical Engineering, University of
Rochester, NY, USA  Faculty of Medicine and Institute of Clinical Radiology,
Ludwig Maximilian University, Munich, Germany  Department of Electrical and
Computer Engineering, University of Rochester, NY, USA  Department of Imaging
Sciences, University of Rochester, NY, USA  Department of Biomedical
Engineering, University of Rochester, NY, USA  Faculty of Medicine and
Institute of Clinical Radiology, Ludwig Maximilian University, Munich, Germany
M. Ali Vosoughi  Department of Electrical and Computer Engineering, University
of Rochester, NY, USA

###### Abstract

We introduce large-scale Augmented Granger Causality (lsAGC) as a method for
connectivity analysis in complex systems. The lsAGC algorithm combines
dimension reduction with source time-series augmentation and uses predictive
time-series modeling for estimating directed causal relationships among time-
series. This method is a multivariate approach, since it is capable of
identifying the influence of each time-series on any other time-series in the
presence of all other time-series of the underlying dynamic system. We
quantitatively evaluate the performance of lsAGC on synthetic directional
time-series networks with known ground truth. As a reference method, we
compare our results with cross-correlation, which is typically used as a
standard measure of connectivity in the functional MRI (fMRI) literature.
Using extensive simulations for a wide range of time-series lengths and two
different signal-to-noise ratios of 5 and 15 dB, lsAGC consistently
outperforms cross-correlation at accurately detecting network connections,
using Receiver Operator Characteristic Curve (ROC) analysis, across all tested
time-series lengths and noise levels. In addition, as an outlook to possible
clinical application, we perform a preliminary qualitative analysis of
connectivity matrices for fMRI data of Autism Spectrum Disorder (ASD) patients
and typical controls, using a subset of 59 subjects of the Autism Brain
Imaging Data Exchange II (ABIDE II) data repository. Our results suggest that
lsAGC, by extracting sparse connectivity matrices, may be useful for network
analysis in complex systems, and may be applicable to clinical fMRI analysis
in future research, such as targeting disease-related classification or
regression tasks on clinical data.

###### keywords:

machine learning, resting-state fMRI, large-scale Augmented Granger Causality,
functional connectivity, autism spectrum disorder

##  1 INTRODUCTION

Currently, the quantification of directed information transfer between
interacting brain areas is one of the most challenging methodological problems
in computational neuroscience. A fundamental problem is identifying
connectivity in very high-dimensional systems. A common practice has been to
transform a high-dimensional system into a simplified representation, e.g. by
clustering, Principal, or Independent Component Analysis. The drawback of such
methodology is that an identified interaction between such simplified
components cannot readily be transferred back into the original high-
dimensional space. Thus, directed interactions between the original network
nodes can no longer be revealed. Although this significantly limits the
interpretation of brain network activities in physiological and disease
states, surprisingly little effort has been devoted to circumvent the
inevitable information loss induced by the aforementioned frequently employed
techniques.

Various methods have been proposed to obtain directional relationships in
multivariate time-series data, e.g., transfer entropy [1] and mutual
information [2]. However, as the multivariate problem’s dimensions increase,
computation of the density function becomes computationally expensive [3, 4].
Under the Gaussian assumption, transfer entropy is equivalent to Granger
causality [5]. However, the computation of multivariate Granger causality for
short time series in large-scale problems is challenging [6, 7]. To address
these problems, we have previously proposed a method for multivariate Granger
causality analysis using linear multivariate auto-regressive (MVAR) modeling,
which simultaneously circumvents the drawbacks of above mentioned
simplification strategies by introducing an invertible dimension reduction
followed by a back-projection of prediction residuals into the original data
space (large-scale Granger Causality, lsGC) [8]. We have also demonstrated the
applicability of this approach to resting-state fMRI analysis [9]. Recently,
we have also presented an alternative multivariate Granger causality analysis
method, large-scale Extended Granger Causality (lsXGC), that uses an augmented
dimension-reduced time-series representation for predicting target time-series
in the original high-dimensional system directly, i.e., without inverting the
dimensionality reduction step [10].

In this paper, we introduce a hybrid of both methods, large-scale Augmented
Granger Causality (lsAGC) that combines both invertible dimension reduction
and time-series augmentation. It first uses an augmented dimension-reduced
time-series representation for prediction in the low-dimensional space,
followed by an inversion of the initial dimension reduction step. In the
following, we explain the lsAGC algorithm and present quantitative results on
synthetic time-series data with known connectivity ground truth. Finally, as
an outlook to possible clinical application, we perform a preliminary
qualitative analysis of connectivity matrices for fMRI data of Autism Spectrum
Disorder (ASD) patients and typical controls, using a subset of the Autism
Brain Imaging Data Exchange II (ABIDE II) data repository.
This work is embedded in our group’s endeavor to expedite artificial
intelligence in biomedical imaging by means of advanced pattern recognition
and machine learning methods for computational radiology and radiomics, e.g.,
[ 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68
].

##  2 DATA

###  2.1 Synthetic Networks with Known Ground Truth for Quantitative Analysis

We quantitatively evaluate the performance of lsAGC for network structure
recovery using synthetic networks with known ground truth. We constructed
ground truth networks with N = 50 nodes, each containing 5 modules of 10 nodes
with high (low) probability for the existence of directed intra- (inter-)
module connections. We simulated two values of additive white Gaussian noise,
with signal-to-noise ratios (SNR) of 15 dB and 5 dB, and repeated the
experiment 100 times with different noise seed. Networks were realized as
noisy stationary multivariate auto-regressive (MVAR) processes of model order
p = 2 in each of T = 1000 temporal process iterations. The network structure
was adapted from [69] and [56].

###  2.2 Functional MRI Data for Qualitative Analysis

The following explanation of participants and data in this section follows the
description in [70]: The Autism Brain Imaging Data Exchange II (ABIDE II)
initiative has made publicly available MRI data from multiple sites. We used
resting-state fMRI data from 59 participants of the online ABIDE II repository
(http:// fcon_1000.projects.nitrc.org/indi/abide) in this analysis, namely the
data from the Olin Neuropsychiatry Research Center, Institute of Living at
Hartford Hospital [71].

This data set consists of 24 ASD subjects (18-31 years) and 35 typical
controls (18-31 years). Autism diagnosis was based on the ASD cutoff of the
Autism Diagnostic Observation Schedule-Generic (ADOS-G). The typical controls
(TC) were screened for autism using the ADOS-G and any psychiatric disorder
based on the Structured Clinical Interview for DSM-IV Axis I Disorders-
Research Version (SCID-I RV) [71].

In this data set, resting-state fMRI scans were obtained from all subjects
using Siemens Magnetom Skyra Syngo MR D13. The study protocol included: (i)
High-resolution structural imaging using T1-weighted magnetization-prepared
rapid gradient-echo (MPRAGE) sequence, TE = 2.88 ms, TR = 2200 ms, isotropic
voxel size 1 mm, flip angle $13^{\circ}$. (ii) Resting-state fMRI scans with
TE = 30 ms, TR = 475 ms, flip angle $60^{\circ}$, isotropic voxel size of 3
mm. The acquisition lasted 7 minutes and 37 seconds.

The fMRI data used in this study were pre-processed using standard
methodology. Motion correction, brain extraction, and correction for slice
timing acquisition were performed. Additional nuisance regression for removing
variations due to head motion and physiological processes was carried out.
Each data set was finally registered to the 2 mm MNI standard space using a
12-parameter affine transformation. The functional data used in this study was
pre-processed using the CONN toolbox (www.nitrc.org/projects/conn) [72].
Additionally, the time-series were normalized to zero mean and unit standard
deviation to focus on signal dynamics rather than amplitude [73]. Finally, the
brain was parcellated into 90 regions as defined by the Automated Anatomic
Labelling (AAL) template [74]. Each regional time-series was represented by
the average time-series of all voxels included in the region.

##  3 ALGORITHM

Large-scale Augmented Granger Causality (lsAGC) has been developed based on 1)
the principle of original Granger causality, which quantifies the causal
influence of time-series $\mathbf{x_{s}}$ on time-series $\mathbf{x_{t}}$ by
quantifying the amount of improvement in the prediction of $\mathbf{x_{t}}$ in
presence of $\mathbf{x_{s}}$. 2) the idea of dimension reduction, which
resolves the problem of the tackling a under-determined system, which is
frequently faced in fMRI analysis, since the number of acquired temporal
samples usually is not sufficient for estimating the model parameters [9].

Consider the ensemble of time-series $\mathcal{X}\in\mathbb{R}^{N\times T}$,
where $N$ is the number of time-series (Regions Of Interest – ROIs) and $T$
the number of temporal samples. Let
$\mathcal{X}=(\mathbf{x_{1}},\mathbf{x_{2}},\dots,\mathbf{x_{N}})^{\mathsf{T}}$
be the whole multidimensional system and $x_{i}\in\mathbb{R}^{1\times T}$ a
single time-series with $i=1,2,\dots,N$, where
$\mathbf{x_{i}}=(x_{i}(1),x_{i}(2),\dots,x_{i}(T))$. In order to overcome the
under-determined problem, first $\mathcal{X}$ will be decomposed into its
first $p$ high-variance principal components
$\mathcal{Z}\in\mathbb{R}^{p\times T}$ using Principal Component Analysis
(PCA), i.e.,

$\mathcal{Z}=W\mathcal{X},$ (1)

where $W\in\mathbb{R}^{p\times N}$ represents the PCA coefficient matrix.
Subsequently, the dimension-reduced time-series ensemble $\mathcal{Z}$ is
augmented by one original time-series $\mathbf{x_{s}}$ yielding a dimension-
reduced augmented time-series ensemble $\mathcal{Y}\in\mathbb{R}^{(p+1)\times
T}$ for estimating the influence of $\mathbf{x_{s}}$ on all other time-series.

Following this, we locally predict the dimension-reduced representation
$\mathcal{Z}$ of the original high-dimensional system $\mathcal{X}$ at each
time sample $t$, i.e. $\mathcal{Z}(t)\in\mathbb{R}^{p\times 1}$ by calculating
an estimate $\hat{\mathcal{Z}}_{\mathbf{x_{s}}}(t)$. To this end, we fit an
affine model based on a vector of $m$ vector of m time samples of
$\mathcal{Y}(\tau)\in\mathbb{R}^{(p+1)\times 1}$($\tau=t-1,t-2,\dots,t-m$),
which is $\mathbf{y}(t)\in\mathbb{R}^{m.(p+1)\times 1}$, and a parameter
matrix $\mathcal{A}\in\mathbb{R}^{p\times m.(p+1)}$ and a constant bias vector
$\mathbf{b}\in\mathbb{R}^{p\times 1}$,

$\hat{\mathcal{Z}}_{\mathbf{x_{s}}}(t)=\mathcal{A}\mathbf{y}(t)+\mathbf{b},~{}~{}t=m+1,m+2,\dots,T.$
(2)

Subsequently, we use the prediction $\hat{\mathcal{Z}}_{\mathbf{x_{s}}}(t)$ to
calculate an estimate of $\mathcal{X}$ at time $t$, i.e.
$\mathcal{X}(t)\in\mathbb{R}^{N\times 1}$ by inverting the PCA of equation
(2), i.e.

$\mathcal{X}=W^{\dagger}\mathcal{Z},$ (3)

where $W^{\dagger}\in\mathbb{R}^{N\times p}$ represents the inverse of the PCA
coefficient matrix $W$, which is calculated as the Moore–Penrose pseudoinverse
of $W$.

Now $\hat{\mathcal{X}}_{\setminus{\mathbf{x_{s}}}}(t)$, which is the
prediction of $\mathcal{X}(t)$ without the information of $\mathbf{x_{s}}$,
will be estimated. The estimation processes is identical to the previous one,
with the only difference being that we have to remove the augmented time-
series $\mathbf{x_{s}}$ and its corresponding column in the PCA coefficient
matrix $W$.

The computation of a lsAGC index is based on comparing the variance of the
prediction errors obtained with and without consideration of $\mathbf{x_{s}}$.
The lsAGC index $f_{\mathbf{x_{s}}\xrightarrow{}\mathbf{x_{t}}}$ , which
indicates the influence of $\mathbf{x_{s}}$ on $\mathbf{x_{t}}$, can be
calculated by the following equation:

$f_{\mathbf{x_{s}}\xrightarrow{}\mathbf{x_{t}}}=\log{\frac{\mathrm{var}(e_{s})}{\mathrm{var}(e_{\setminus
s})}},$ (4)

where $e_{\setminus s}$ is the error in predicting $\mathbf{x_{t}}$ when
$\mathbf{x_{s}}$ was not considered, and $e_{s}$ is the error, when
$\mathbf{x_{s}}$ was used. Based on preliminary analyses, in this study, we
set $p=7$ and $m=3$.

##  4 RESULTS

Quantitative Analysis of Synthetic Networks with Known Ground Truth: Network
reconstruction results for the synthetic networks with known ground truth,
using the Area Under the Curve (AUC) for Receiver Operating Characteristic
(ROC) analysis, are shown in Fig. 1. For each time-series length and each
noise level, we performed 100 simulations. As can be seen from Fig. 1, lsAGC
consistently outperforms cross-correlation in its ability to accurately
recover network structure over a wide range of time-series-lengths in both
high- and low-noise scenarios, with a mean AUC for lsAGC for a time-series
length of 1000 temporal samples equal to 98.9% and 97.1% for signal-to-noise
values of 15 and 5 dB, respectively. On the other hand, cross-correlation
performs quite poorly compared to lsAGC with its mean AUC ranging around 0.5
for all examined time-series lengths and noise levels, equivalent to the
quality of randomly guessing the presence or absence of network connections.

Figure 1: Quantitative performance comparison of cross-correlation and lsAGC
for recovery of synthetic networks. The vertical axis is the Area Under the
Curve (AUC) for Receiver Operating Characteristics (ROC) analysis, where an
AUC = 1 indicates a perfect network recovery and AUC = 0.5 random assignment.
Whiskers are related to the 95% confidence interval, green diamonds represent
outliers, orange lines represent medians, and boxes are drawn from the first
quartile to the third quartile. It is clearly seen that lsAGC outperforms
cross-correlation over all tested time-series lengths and noise levels.

Qualitative Analysis of Connectivity Matrices Extracted from fMRI Data:
Averaged connectivity matrices, which were extracted using lsAGC and cross-
correlation, are shown in Fig. 2 for both healthy controls and ASD patients.
These matrices were obtained by calculating and then averaging over the
connectivity matrices of the 24 ASD patients and 35 typical controls, using
the proposed lsAGC algorithm as well as conventional cross-correlation
analysis. Visual inspection of the mean connectivity matrices in Fig. 2
reveals subtle differences between ASD patients and healthy controls for both
methods, which may be exploited for classification among the two cohorts in
future research. We also find from visual inspection of Fig. 2 that the
features extracted by the two methods are likely different, where the mean
connectivity matrices for lsAGC appear to be more “sparse” than for cross-
correlation. To quantify this qualitative visual impression, we calculated the
entropy of the connectivity matrix elements for each of the 59 subjects as a
surrogate for matrix “sparseness”. The mean entropy for healthy controls with
lsAGC and correlation was 1.66 $\pm$ 0.14 and 4.76 $\pm$ 0.07, respectively,
and the mean entropy for the ASD patients with lsAGC and correlation was 1.75
$\pm$ 0.12 and 4.78 $\pm$ 0.07, respectively. I.e., for both cohorts, the
sparseness of lsAGC connectivity matrices appears to be higher than for cross-
correlation analysis. We found that this difference between methods, as
expressed by the entropy of the connectivity matrix elements, was
statistically significant (Mann-Whitney U-test, $p<10^{-8}$). We conclude that
lsAGC may be useful for disease-related classification or regression tasks on
clinical fMRI data, because it may extract relevant features potentially not
captured by cross-correlation, which is currently used as the mainstay of fMRI
connectivity analysis. This hypothesis can be further investigated in future
research.

Figure 2: Averaged connectivity matrices: top left: average connectivity
matrix of healthy control subjects using lsAGC, top right: average
connectivity matrix of ASD patients using lsAGC, bottom left: average
connectivity matrix of healthy control subjects using cross-correlation, and
bottom right: average connectivity matrix of ASD patients using cross-
correlation. Note that the different methods capture different connectivity
features, and that there are slight differences of connectivity patterns
between healthy subjects and ASD patients. Also, the lsAGC connectivity
matrices appear to be significantly “sparser” than cross-correlation matrices.
This observation is quantitatively confirmed by calculating the entropy over
the matrix elements, as explained in the text.

##  5 CONCLUSIONS

In this work, we have introduced large-scale Augmented Granger Causality
(lsAGC) as a method for connectivity analysis in complex systems. The lsAGC
algorithm combines dimension reduction with source time-series augmentation
and uses multivariate predictive time-series modeling for estimating directed
causal relationships among time-series. We quantitatively evaluated the
performance of lsAGC on synthetic directional time-series networks with known
ground truth. Using simulations for a wide range of time-series lengths and
different signal-to-noise ratios, we compared lsAGC with cross-correlation,
which is currently used as the clinical standard for fMRI connectivity
analysis. We found that lsAGC consistently outperformed cross-correlation at
accurately detecting network connections. In addition, we performed a
preliminary qualitative analysis of connectivity matrices for fMRI data of
Autism Spectrum Disorder (ASD) patients and typical controls, using a subset
of the ABIDE II data repository. Our results suggest that lsAGC, by extracting
sparse connectivity matrices, may be useful for network analysis in complex
systems, and may be applicable to clinical fMRI analysis in future research,
such as targeting disease-related classification or regression tasks on
clinical data.

###### Acknowledgements.

This research was funded by Ernest J. Del Monte Institute for Neuroscience
Award from the Harry T. Mangurian Jr. Foundation. This work was conducted as a
Practice Quality Improvement (PQI) project related to American Board of
Radiology (ABR) Maintenance of Certificate (MOC) for Prof. Dr. Axel Wismüller.
This work is not being and has not been submitted for publication or
presentation elsewhere.

## References

  * [1] Schreiber, T., “Measuring information transfer,” Physical review letters 85(2), 461 (2000).
  * [2] Kraskov, A., Stögbauer, H., and Grassberger, P., “Estimating mutual information,” Physical review E 69(6), 066138 (2004).
  * [3] Mozaffari, M. and Yilmaz, Y., “Online multivariate anomaly detection and localization for high-dimensional settings,” arXiv preprint arXiv:1905.07107 (2019).
  * [4] Mozaffari, M. and Yilmaz, Y., “Online anomaly detection in multivariate settings,” in [2019 IEEE 29th International Workshop on Machine Learning for Signal Processing (MLSP) ], 1–6, IEEE (2019).
  * [5] Barnett, L., Barrett, A. B., and Seth, A. K., “Granger causality and transfer entropy are equivalent for Gaussian variables,” Physical review letters 103(23), 238701 (2009).
  * [6] Vosoughi, M. A. and Wismüller, A., “Large-scale kernelized Granger causality to infer topology of directed graphs with applications to brain networks,” arXiv preprint arXiv:2011.08261 (2020).
  * [7] Wismüller, A., DSouza, A. M., Abidin, A. Z., and Vosoughi, M. A., “Large-scale nonlinear Granger causality: A data-driven, multivariate approach to recovering directed networks from short time-series data,” arXiv preprint arXiv:2009.04681 (2020).
  * [8] D’Souza, A. M., Abidin, A. Z., Leistritz, L., and Wismüller, A., “Large-scale Granger causality analysis on resting-state functional MRI,” in [Medical Imaging 2016: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 9788, 97880L, International Society for Optics and Photonics (2016).
  * [9] DSouza, A. M., Abidin, A. Z., Leistritz, L., and Wismüller, A., “Exploring connectivity with large-scale Granger causality on resting-state functional MRI,” Journal of neuroscience methods 287, 68–79 (2017).
  * [10] Vosoughi, M. A. and Wismüller, A., “Large-scale extended Granger causality for classification of marijuana users from functional MRI,” arXiv preprint arXiv:2101.01832 (2021).
  * [11] Leinsinger, G., Schlossbauer, T., Scherr, M., Lange, O., Reiser, M., and Wismüller, A., “Cluster analysis of signal-intensity time course in dynamic breast MRI: does unsupervised vector quantization help to evaluate small mammographic lesions?,” European radiology 16(5), 1138–1146 (2006).
  * [12] Wismüller, A., Vietze, F., Behrends, J., Meyer-Baese, A., Reiser, M., and Ritter, H., “Fully automated biomedical image segmentation by self-organized model adaptation,” Neural Networks 17(8-9), 1327–1344 (2004).
  * [13] Nattkemper, T. W. and Wismüller, A., “Tumor feature visualization with unsupervised learning,” Medical Image Analysis 9(4), 344–351 (2005).
  * [14] Bunte, K., Hammer, B., Wismüller, A., and Biehl, M., “Adaptive local dissimilarity measures for discriminative dimension reduction of labeled data,” Neurocomputing 73(7-9), 1074–1092 (2010).
  * [15] Wismüller, A., Vietze, F., and Dersch, D. R., “Segmentation with neural networks,” in [Handbook of medical imaging ], 107–126, Academic Press, Inc. (2000).
  * [16] Hoole, P., Wismüller, A., Leinsinger, G., Kroos, C., Geumann, A., and Inoue, M., “Analysis of tongue configuration in multi-speaker, multi-volume MRI data,” (2000).
  * [17] Wismüller, A., “Exploratory morphogenesis (XOM): a novel computational framework for self-organization,” Ph. D. thesis, Technical University of Munich, Department of Electrical and Computer Engineering (2006).
  * [18] Wismüller, A., Dersch, D. R., Lipinski, B., Hahn, K., and Auer, D., “A neural network approach to functional MRI pattern analysis—clustering of time-series by hierarchical vector quantization,” in [International Conference on Artificial Neural Networks ], 857–862, Springer (1998).
  * [19] Wismüller, A., Vietze, F., Dersch, D. R., Behrends, J., Hahn, K., and Ritter, H., “The deformable feature map-a novel neurocomputing algorithm for adaptive plasticity in pattern analysis,” Neurocomputing 48(1-4), 107–139 (2002).
  * [20] Behrends, J., Hoole, P., Leinsinger, G. L., Tillmann, H. G., Hahn, K., Reiser, M., and Wismüller, A., “A segmentation and analysis method for MRI data of the human vocal tract,” in [Bildverarbeitung für die Medizin 2003 ], 186–190, Springer (2003).
  * [21] Wismüller, A., “Neural network computation in biomedical research: chances for conceptual cross-fertilization,” Theory in Biosciences (1997).
  * [22] Bunte, K., Hammer, B., Villmann, T., Biehl, M., and Wismüller, A., “Exploratory observation machine (XOM) with Kullback-Leibler divergence for dimensionality reduction and visualization.,” in [ESANN ], 10, 87–92 (2010).
  * [23] Wismüller, A., Vietze, F., Dersch, D. R., Hahn, K., and Ritter, H., “The deformable feature map—adaptive plasticity for function approximation,” in [International Conference on Artificial Neural Networks ], 123–128, Springer (1998).
  * [24] Wismüller, A., “The exploration machine–a novel method for data visualization,” in [International Workshop on Self-Organizing Maps ], 344–352, Springer (2009).
  * [25] Wismüller, A., “Method, data processing device and computer program product for processing data,” (July 28 2009).  US Patent 7,567,889.
  * [26] Huber, M. B., Nagarajan, M., Leinsinger, G., Ray, L. A., and Wismüller, A., “Classification of interstitial lung disease patterns with topological texture features,” in [Medical Imaging 2010: Computer-Aided Diagnosis ], 7624, 762410, International Society for Optics and Photonics (2010).
  * [27] Wismüller, A., “The exploration machine: a novel method for analyzing high-dimensional data in computer-aided diagnosis,” in [Medical Imaging 2009: Computer-Aided Diagnosis ], 7260, 72600G, International Society for Optics and Photonics (2009).
  * [28] Bunte, K., Hammer, B., Villmann, T., Biehl, M., and Wismüller, A., “Neighbor embedding XOM for dimension reduction and visualization,” Neurocomputing 74(9), 1340–1350 (2011).
  * [29] Meyer-Bäse, A., Lange, O., Wismüller, A., and Ritter, H., “Model-free functional MRI analysis using topographic independent component analysis,” International journal of neural systems 14(04), 217–228 (2004).
  * [30] Wismüller, A., “A computational framework for nonlinear dimensionality reduction and clustering,” in [International Workshop on Self-Organizing Maps ], 334–343, Springer (2009).
  * [31] Meyer-Base, A., Auer, D., and Wismüller, A., “Topographic independent component analysis for fMRI signal detection,” in [Proceedings of the International Joint Conference on Neural Networks, 2003. ], 1, 601–605, IEEE (2003).
  * [32] Meyer-Baese, A., Schlossbauer, T., Lange, O., and Wismüller, A., “Small lesions evaluation based on unsupervised cluster analysis of signal-intensity time courses in dynamic breast MRI,” International journal of biomedical imaging 2009 (2009).
  * [33] Wismüller, A., Lange, O., Auer, D., and Leinsinger, G., “Model-free functional MRI analysis for detecting low-frequency functional connectivity in the human brain,” in [Medical Imaging 2010: Computer-Aided Diagnosis ], 7624, 76241M, International Society for Optics and Photonics (2010).
  * [34] Meyer-Bäse, A., Saalbach, A., Lange, O., and Wismüller, A., “Unsupervised clustering of fMRI and MRI time series,” Biomedical Signal Processing and Control 2(4), 295–310 (2007).
  * [35] Huber, M. B., Nagarajan, M. B., Leinsinger, G., Eibel, R., Ray, L. A., and Wismüller, A., “Performance of topological texture features to classify fibrotic interstitial lung disease patterns,” Medical Physics 38(4), 2035–2044 (2011).
  * [36] Wismüller, A., Verleysen, M., Aupetit, M., and Lee, J. A., “Recent advances in nonlinear dimensionality reduction, manifold and topological learning.,” in [ESANN ], (2010).
  * [37] Meyer-Baese, A., Lange, O., Wismüller, A., and Hurdal, M. K., “Analysis of dynamic susceptibility contrast MRI time series based on unsupervised clustering methods,” IEEE Transactions on Information Technology in Biomedicine 11(5), 563–573 (2007).
  * [38] Wismüller, A., Behrends, J., Hoole, P., Leinsinger, G. L., Reiser, M. F., and Westesson, P.-L., “Human vocal tract analysis by in vivo 3d MRI during phonation: a complete system for imaging, quantitative modeling, and speech synthesis,” in [International Conference on Medical Image Computing and Computer-Assisted Intervention ], 306–312, Springer (2008).
  * [39] Wismüller, A., “Method and device for representing multichannel image data,” (Nov. 17 2015).  US Patent 9,189,846.
  * [40] Huber, M. B., Bunte, K., Nagarajan, M. B., Biehl, M., Ray, L. A., and Wismüller, A., “Texture feature ranking with relevance learning to classify interstitial lung disease patterns,” Artificial intelligence in medicine 56(2), 91–97 (2012).
  * [41] Wismüller, A., Meyer-Baese, A., Lange, O., Reiser, M. F., and Leinsinger, G., “Cluster analysis of dynamic cerebral contrast-enhanced perfusion MRI time-series,” IEEE transactions on medical imaging 25(1), 62–73 (2005).
  * [42] Twellmann, T., Saalbach, A., Muller, C., Nattkemper, T. W., and Wismüller, A., “Detection of suspicious lesions in dynamic contrast enhanced MRI data,” in [The 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society ], 1, 454–457, IEEE (2004).
  * [43] Otto, T. D., Meyer-Baese, A., Hurdal, M., Sumners, D., Auer, D., and Wismüller, A., “Model-free functional MRI analysis using cluster-based methods,” in [Intelligent Computing: Theory and Applications ], 5103, 17–24, International Society for Optics and Photonics (2003).
  * [44] Varini, C., Nattkemper, T. W., Degenhard, A., and Wismüller, A., “Breast MRI data analysis by lle,” in [2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541) ], 3, 2449–2454, IEEE (2004).
  * [45] Huber, M. B., Lancianese, S. L., Nagarajan, M. B., Ikpot, I. Z., Lerner, A. L., and Wismüller, A., “Prediction of biomechanical properties of trabecular bone in mr images with geometric features and support vector regression,” IEEE Transactions on Biomedical Engineering 58(6), 1820–1826 (2011).
  * [46] Meyer-Base, A., Pilyugin, S. S., and Wismüller, A., “Stability analysis of a self-organizing neural network with feedforward and feedback dynamics,” in [2004 IEEE International Joint Conference on Neural Networks (IEEE Cat. No. 04CH37541) ], 2, 1505–1509, IEEE (2004).
  * [47] Meyer-Baese, A., Lange, O., Schlossbauer, T., and Wismüller, A., “Computer-aided diagnosis and visualization based on clustering and independent component analysis for breast MRI,” in [2008 15th IEEE International Conference on Image Processing ], 3000–3003, IEEE (2008).
  * [48] Wismüller, A., Meyer-Bäse, A., Lange, O., Schlossbauer, T., Kallergi, M., Reiser, M., and Leinsinger, G., “Segmentation and classification of dynamic breast magnetic resonance image data,” Journal of Electronic Imaging 15(1), 013020 (2006).
  * [49] Bhole, C., Pal, C., Rim, D., and Wismüller, A., “3d segmentation of abdominal ct imagery with graphical models, conditional random fields and learning,” Machine vision and applications 25(2), 301–325 (2014).
  * [50] Nagarajan, M. B., Coan, P., Huber, M. B., Diemoz, P. C., Glaser, C., and Wismüller, A., “Computer-aided diagnosis in phase contrast imaging x-ray computed tomography for quantitative characterization of ex vivo human patellar cartilage,” IEEE Transactions on Biomedical Engineering 60(10), 2896–2903 (2013).
  * [51] Wismüller, A., Meyer-Bäse, A., Lange, O., Auer, D., Reiser, M. F., and Sumners, D., “Model-free functional MRI analysis based on unsupervised clustering,” Journal of Biomedical Informatics 37(1), 10–18 (2004).
  * [52] Meyer-Baese, A., Wismüller, A., Lange, O., and Leinsinger, G., “Computer-aided diagnosis in breast MRI based on unsupervised clustering techniques,” in [Intelligent Computing: Theory and Applications II ], 5421, 29–37, International Society for Optics and Photonics (2004).
  * [53] Nagarajan, M. B., Coan, P., Huber, M. B., Diemoz, P. C., Glaser, C., and Wismüller, A., “Computer-aided diagnosis for phase-contrast x-ray computed tomography: quantitative characterization of human patellar cartilage with high-dimensional geometric features,” Journal of digital imaging 27(1), 98–107 (2014).
  * [54] Nagarajan, M. B., Huber, M. B., Schlossbauer, T., Leinsinger, G., Krol, A., and Wismüller, A., “Classification of small lesions on dynamic breast MRI: Integrating dimension reduction and out-of-sample extension into cadx methodology,” Artificial intelligence in medicine 60(1), 65–77 (2014).
  * [55] Yang, C.-C., Nagarajan, M. B., Huber, M. B., Carballido-Gamio, J., Bauer, J. S., Baum, T. H., Eckstein, F., Lochmüller, E.-M., Majumdar, S., Link, T. M., et al., “Improving bone strength prediction in human proximal femur specimens through geometrical characterization of trabecular bone microarchitecture and support vector regression,” Journal of electronic imaging 23(1), 013013 (2014).
  * [56] Wismüller, A., Nagarajan, M. B., Witte, H., Pester, B., and Leistritz, L., “Pair-wise clustering of large scale Granger causality index matrices for revealing communities,” in [Medical Imaging 2014: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 9038, 90381R, International Society for Optics and Photonics (2014).
  * [57] Wismüller, A., Wang, X., DSouza, A. M., and Nagarajan, M. B., “A framework for exploring non-linear functional connectivity and causality in the human brain: mutual connectivity analysis (mca) of resting-state functional MRI with convergent cross-mapping and non-metric clustering,” arXiv preprint arXiv:1407.3809 (2014).
  * [58] Schmidt, C., Pester, B., Nagarajan, M., Witte, H., Leistritz, L., and Wismüller, A., “Impact of multivariate Granger causality analyses with embedded dimension reduction on network modules,” in [2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society ], 2797–2800, IEEE (2014).
  * [59] Wismüller, A., Abidin, A. Z., D’Souza, A. M., Wang, X., Hobbs, S. K., Leistritz, L., and Nagarajan, M. B., “Nonlinear functional connectivity network recovery in the human brain with mutual connectivity analysis (MCA): convergent cross-mapping and non-metric clustering,” in [Medical Imaging 2015: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 9417, 94170M, International Society for Optics and Photonics (2015).
  * [60] Wismüller, A., Abidin, A. Z., DSouza, A. M., and Nagarajan, M. B., “Mutual connectivity analysis (MCA) for nonlinear functional connectivity network recovery in the human brain using convergent cross-mapping and non-metric clustering,” in [Advances in Self-Organizing Maps and Learning Vector Quantization ], 217–226, Springer (2016).
  * [61] Schmidt, C., Pester, B., Schmid-Hertel, N., Witte, H., Wismüller, A., and Leistritz, L., “A multivariate Granger causality concept towards full brain functional connectivity,” PloS one 11(4) (2016).
  * [62] Abidin, A. Z., Chockanathan, U., DSouza, A. M., Inglese, M., and Wismüller, A., “Using large-scale Granger causality to study changes in brain network properties in the clinically isolated syndrome (CIS) stage of multiple sclerosis,” in [Medical Imaging 2017: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 10137, 101371B, International Society for Optics and Photonics (2017).
  * [63] Chen, L., Wu, Y., DSouza, A. M., Abidin, A. Z., Wismüller, A., and Xu, C., “MRI tumor segmentation with densely connected 3d cnn,” in [Medical Imaging 2018: Image Processing ], 10574, 105741F, International Society for Optics and Photonics (2018).
  * [64] Abidin, A. Z., DSouza, A. M., Nagarajan, M. B., Wang, L., Qiu, X., Schifitto, G., and Wismüller, A., “Alteration of brain network topology in HIV-associated neurocognitive disorder: A novel functional connectivity perspective,” NeuroImage: Clinical 17, 768–777 (2018).
  * [65] Abidin, A. Z., Deng, B., DSouza, A. M., Nagarajan, M. B., Coan, P., and Wismüller, A., “Deep transfer learning for characterizing chondrocyte patterns in phase contrast x-ray computed tomography images of the human patellar cartilage,” Computers in biology and medicine 95, 24–33 (2018).
  * [66] DSouza, A. M., Abidin, A. Z., Chockanathan, U., Schifitto, G., and Wismüller, A., “Mutual connectivity analysis of resting-state functional MRI data with local models,” NeuroImage 178, 210–223 (2018).
  * [67] Chockanathan, U., DSouza, A. M., Abidin, A. Z., Schifitto, G., and Wismüller, A., “Automated diagnosis of HIV-associated neurocognitive disorders using large-scale Granger causality analysis of resting-state functional MRI,” Computers in Biology and Medicine 106, 24–30 (2019).
  * [68] Wismüller, A. and Vosoughi, M. A., “Classification of schizophrenia from functional MRI using large-scale extended Granger causality,” arXiv preprint (2021).
  * [69] Baccalá, L. A. and Sameshima, K., “Partial directed coherence: a new concept in neural structure determination,” Biological cybernetics 84(6), 463–474 (2001).
  * [70] DSouza, A. M., Abidin, A. Z., and Wismüller, A., “Classification of autism spectrum disorder from resting-state fMRI with mutual connectivity analysis,” in [Medical Imaging 2019: Biomedical Applications in Molecular, Structural, and Functional Imaging ], 10953, 109531D, International Society for Optics and Photonics (2019).
  * [71] First, M., Spitzer, R., Gibbon, M., and Williams, J., “Structured clinical interview for dsm-iv axis i disorders (scid). new york, new york state psychiatric institute,” Biometrics Research (1995).
  * [72] Whitfield-Gabrieli, S. and Nieto-Castanon, A., “Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks,” Brain connectivity 2(3), 125–141 (2012).
  * [73] Wismüller, A., Lange, O., Dersch, D. R., Leinsinger, G. L., Hahn, K., Pütz, B., and Auer, D., “Cluster analysis of biomedical image time-series,” International Journal of Computer Vision 46(2), 103–128 (2002).
  * [74] Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B., and Joliot, M., “Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni MRI single-subject brain,” Neuroimage 15(1), 273–289 (2002).


# Exact exponential synchronization rate of high-dimensional Kuramoto models
with identical oscillators and digraphs

Shanshan Peng<EMAIL_ADDRESS>Jinxing Zhang<EMAIL_ADDRESS>Jiandong Zhu<EMAIL_ADDRESS>Jianquan Lu<EMAIL_ADDRESS>Xiaodi Li
<EMAIL_ADDRESS>School of Mathematical Sciences, Nanjing Normal University,
Nanjing, 210023, PRC  School of Mathematics, Southeast University, Nanjing,
211189, PRC  School of Mathematics and Statistics, Shandong Normal University,
Jinan, 250358, PRC

###### Abstract

For the high-dimensional Kuramoto model with identical oscillators under a
general digraph that has a directed spanning tree, although exponential
synchronization was proved under some initial state constraints, the exact
exponential synchronization rate has not been revealed until now. In this
paper, the exponential synchronization rate is precisely determined as the
smallest non-zero real part of Laplacian eigenvalues of the digraph. Our
obtained result extends the existing results from the special case of strongly
connected balanced digraphs to the condition of general digraphs owning
directed spanning trees, which is the weakest condition for synchronization
from the aspect of network structure. Moreover, our adopted method is
completely different from and much more elementary than the previous
differential geometry method.

###### keywords:

High-dimensional Kuramoto model; Exponential synchronization rate; Directed
graph.

††thanks: This work is supported in part by National Natural Science
Foundation (NNSF) of China under Grants 61673012 and 11971240.

, , , ,

##  1 Introduction

The well-known Kuramoto model proposed by Yoshiki Kuramoto in [8] is one of
the most successful mathematical models to describe collective behaviors of
complex dynamical networks. Kuramoto model and its various generalized forms
have many applications in engineering, neuroscience, physics and so on [18].
For Kuramoto model, exponential synchronization is an typical collective
behavior and an interesting theoretical issue [6], [7], [20], [21]. The
dynamics of the general Kuramoto model is described by

$\begin{split}\dot{\theta}_{i}&=\omega_{i}+k\sum\limits_{j=1}^{m}a_{ij}\sin(\theta_{j}-\theta_{i}),\
\ i=1,2,\dots,m,\end{split}$ (1)

where the $m\times m$ matrix $(a_{ij})$ is the adjacency matrix of the
interconnecting graph, $\theta_{i}$ is the $i$-th oscillator’s phase angle and
$\omega_{i}$ is the natural frequency. The $n$-dimensional vector-valued
Kuramoto model is described as follows:

$\dot{r}_{i}=\Omega_{i}r_{i}+k\sum\limits_{j=1}^{m}a_{ij}(r_{j}-\frac{r_{i}^{T}r_{j}}{r_{i}^{T}r_{i}}r_{i}),\
i=1,2,\dots,m,$ (2)

where $r_{i}\in{\mathbf{R}^{n}}$ is the $i$-th oscillator’s state,
$\Omega_{i}$ is an $n\times n$ skew-symmetric matrix. When $n=2$ and
$r_{i}=[\cos\theta_{i},\ \sin\theta_{i}]^{\mathrm{T}}$, by the $n$-dimensional
Kuramoto model (2), it is easy to deduce the original Kuramoto model (1). If
(2) has the identical oscillators, i.e. $\Omega_{i}=\Omega$ for
$i=1,2,\dots,m$, one can assume, without loss of generality, that
$\Omega_{i}=0$ for $i=1,2,\dots,m$ [24]. Limited on the unit sphere
$\mathbf{S}^{n-1}$, dynamical network (2) with $\Omega_{i}=0$ and $k=1$ is
reduced to

$\dot{r}_{i}=\sum\limits_{j=1}^{m}a_{ij}\big{(}r_{j}-(r_{i}^{T}r_{j})r_{i}\big{)},\
i=1,2,\dots,m,$ (3)

which is first proposed in [16] as a swarm model on spheres, and can be used
to solve the max-cut problem. Due to Lohe’s pioneering literatures such as
[11] and [12], a class of high-dimensional Kuramoto models is also called Lohe
model, which has potential applications to quantum systems. More general
Kuramoto models defined on some matrix manifolds can be seen in [1], [15],
[19].

For the case of complete graphs, many theoretical results on the
synchronization of high-dimensional Kuramoto model have been achieved such as
[5] , [16] and [10]. In our early paper [24], the synchronization limited on
an open half-sphere is proved for the case of general undirected connected
graphs. If the states of particles are not limited on an open half-sphere,
some almost global synchronization results can also be obtained for the case
of undirected connected graphs [14], [25].

For general directed graphs, the synchronization of the high-dimensional
Kuramoto model is proved in [9] by using a differential geometry method.
However, in our paper [23], a simpler approach based on LaSalle invariance
principle is adopted to demonstrate the synchronization under the general
directed graph condition. Just like the original Kuramoto model, high-
dimensional Kuramoto model also has the dynamical property of exponential
synchronization (see [4], [5], [9] and [22]). But so far, the exact
exponential synchronization rate is only obtained for a very special kind of
digraphs, i.e. strongly connected balanced digraphs. In our recent paper [22],
the exponential synchronization is achieved for a general digraph admitting a
spanning tree, but the exact exponential synchronization rate has not yet been
obtained.

In this paper, for the high-dimensional Kuramoto model under a general digraph
containing a spanning tree, it is proved the exponential synchronization rate
is also exactly $\mathrm{Re}(\lambda_{2})$ just like the case of strongly
connected balanced digraphs considered in [9], where $\lambda_{2}$ is the
eigenvalue of the Laplacian $L$ with the smallest non-zero real part.

The rest of this paper is organized as follows. Section 2 gives some
preliminaries and the problem statement. Section 3 includes our main results.
Section 4 is a summary.

##  2 Preliminaries and Problem Statement

For high-dimensional Kuramoto model (3), exponential synchronization is
defined as follows:

######  Definition 2.1.

It is said that the exponential synchronization for (3) is achieved if there
exists a $\mu>0$ and a constant $c>0$ such that

$||r_{i}(t)-r_{j}(t)||\leq c{\rm e}^{-\mu t}\ \ \forall\ i,\ j=1,2,\dots,m.$
(4)

Moreover, we say that the exponential synchronization rate is at least $\mu$
when (4) holds. The minimum $\mu$ satisfying (4) is just called the
exponential synchronization rate.

The first result on the exponential synchronization of $(\ref{eq3})$ under
digraphs is obtained in [9], which is based on differential geometry method.
The main contribution of Theorem 1 and Corollary 1 of [9] can be rewritten as
follows:

######  Proposition 2.1.

Consider the high-dimensional Kuramoto model $(\ref{eq3})$. If the
interconnecting graph is strongly connected and balanced, then the local
exponential synchronization is achieved and the exact exponential
synchronization rate is the smallest non-zero real part of the Laplacian
eigenvalues of the interconnecting digraph.

The target of this paper is to get the exact exponential synchronization rate
under the more general digraph condition and the framework of matrix Riccati
differential equation proposed in our early paper [22].

Let

$e_{ij}=1-r_{i}^{T}r_{j}=\frac{1}{2}||r_{i}-r_{j}||^{2}.$ (5)

It is easily seen that $e_{ij}=e_{ji}$, $e_{ii}=0$ and $0\leq{e_{ij}}\leq{2}$
for every $i,j=1,2,\dots,m$. Let
$E(t)\\!=\\!(e_{ij}(t))\\!\in{\\!\mathbf{R}^{m\times{m}}}$. By [22], the
dynamics of $E(t)$ is described by matrix Riccati differential equation

$\displaystyle\dot{E}$ $\displaystyle\\!=$
$\displaystyle\\!-LE\\!-\\!EL^{\\!\mathrm{T}}\\!\\!-\\!\alpha(E)\mathbf{1}^{\\!\mathrm{T}}\\!\\!-\\!\mathbf{1}\alpha^{\\!\mathrm{T}}\\!(E)\\!+\\!\Lambda(E)\\!E\\!+\\!E\Lambda(E),$
(6)

where $L$ is the Laplacian matrix of the digraph,

$\alpha(E)\\!=\\!(\alpha_{1}(E),\alpha_{2}(E),\dots,\alpha_{m}(E))^{\mathrm{T}}\\!\in{\\!\mathbf{R}^{m}},$
(7)

with each $\alpha_{i}(E)\\!=\\!\sum\limits_{l=1}^{m}a_{il}e_{il}$, and
$\Lambda(E)\\!=\\!\mathrm{diag(}\alpha(E))$ is the diagonal matrix with the
diagonal elements composed of $\alpha_{1}(E)$, $\alpha_{2}(E)$, $\dots$,
$\alpha_{m}(E)$.

For the linear space $\mathbf{R}^{m\times m}$, denote by $S_{m}^{0}$ the
subspace composed of all the symmetric real matrices with all the diagonal
entries being $0$. Then the synchronization error equation (6) is the dynamics
restricted on $S_{m}^{0}$. Let $S_{m}$ and $K_{m}$ be the $m$-order symmetric
matrix subspace and the $m$-order skew-symmetric matrix subspace,
respectively. Then

$\mathbf{R}^{m}=S_{m}\oplus K_{m},\ \ S_{m}=S_{m}^{0}\oplus D_{m},$ (8)

where $D_{m}$ is the $m$-order diagonal matrix subspace and $\oplus$ denotes
the direct sum.

##  3 Main Results

In our early paper [22], Lyapunov functions are designed by using the left
eigenvector of the Laplacian matrix of the digraph and consequently the
exponential synchronization is proved. However, we cannot get the exact
exponential synchronization rate by the approach of Lyapunov functions.

In this section, we turn to the Lyapunov’s first method, i.e. the approximate
exponential linearization method to get the exact synchronization rate. It is
straightforward to check that the approximate linearized system of (6) can be
rewritten as

$\displaystyle\dot{E}$ $\displaystyle\\!=$
$\displaystyle\\!-\big{(}\\!LE\\!+\\!EL^{\\!\mathrm{T}}\\!-\hat{L}\mathrm{vec}(E)\mathbf{1}_{m}^{\\!\mathrm{T}}\\!-\\!\mathbf{1}_{m}(\mathrm{vec}(E))^{\\!\mathrm{T}}\hat{L}^{\\!\mathrm{T}}\\!\big{)},$
(9)

where $\mathrm{vec}(\cdot):\mathbf{R}^{m\times
m}\rightarrow\mathbf{R}^{m^{2}}$ denotes the column-stacking operator,

$\hat{L}\\!=\\!\left[\\!\begin{array}[]{cccc}L_{1}&\\\ &L_{2}\\\ &&\ddots\\\
&&&L_{m}\end{array}\\!\right]\\!,\ \ L_{i}\\!=\\!\mathrm{Row}_{i}(L),\ \
i\\!=\\!1,2,\dots,m.$

Considering (9), we define a linear transformation on $\mathbf{R}^{m\times m}$
by

$\mathcal{T}(X)=LX\\!+\\!XL^{\\!\mathrm{T}}\\!-\hat{L}\mathrm{vec}(X)\mathbf{1}_{m}^{\\!\mathrm{T}}\\!-\\!\mathbf{1}_{m}(\mathrm{vec}(X))^{\\!\mathrm{T}}\hat{L}^{\\!\mathrm{T}}.$
(10)

So the exponential decay rate of $E(t)$ is just the smallest non-zero real
part of eigenvalues of $\mathcal{T}(X)$ restricted on $S_{m}^{0}$. It is not
easy to directly compute the eigenvalues of $\mathcal{T}(X)$ restricted on
$S_{m}^{0}$. We first investigate properties of $\mathcal{T}(X)$ on
$\mathbf{R}^{m}$.

######  Proposition 3.1.

Consider the linear transformation $\mathcal{T}(X)$ defined by (10) and the
Lyapunov mapping $\mathcal{S}(X)=LX+XL^{\mathrm{T}}$. For the subspaces
$S_{m}^{0}$, $D_{m}$ and $K_{m}$, the following statements hold:

(i) if $X\in S_{m}$, then $\mathcal{T}(X)\in S_{m}^{0}$;

(ii) if $X\in K_{m}$, then the projection of $\mathcal{T}(X)$ onto $K_{m}$ is
just $\mathcal{S}(X)$.

(i) If $X=X^{\mathrm{T}}$, it is easy to check that
$\big{(}\mathcal{T}(X)\big{)}^{\mathrm{T}}=\mathcal{T}(X)$. Moreover, a
straightforward computation shows that

$\displaystyle\delta_{i}^{\mathrm{T}}\mathcal{T}(X)\delta_{i}$
$\displaystyle=$ $\displaystyle
2\big{(}\delta_{i}^{\mathrm{T}}LX\delta_{i}-\delta_{i}^{\mathrm{T}}\hat{L}\mathrm{vec}(X)\mathbf{1}_{m}^{\mathrm{T}}\delta_{i}\big{)}$
$\displaystyle=$ $\displaystyle
2\big{(}L_{i}X_{i}-\delta_{i}^{\mathrm{T}}\hat{L}\mathrm{vec}(X)\big{)}=0,$

where $\delta_{i}$ is the $m$-dimensional vector whose $i$-th entry is one and
the others are zero. So $\mathcal{T}(X)\in S_{m}^{0}$.

(ii) If $X=-X^{\mathrm{T}}$, then

$(\mathcal{S}(X))^{\\!\mathrm{T}}=X^{\\!\mathrm{T}}L^{\\!\mathrm{T}}+LX^{\\!\mathrm{T}}=-\mathcal{S}(X).$

Thus $\mathcal{S}(X)\in K_{m}$. From (10), it follows that

$\mathcal{T}(X)=\mathcal{S}(X)-(\hat{L}\mathrm{vec}(X)\mathbf{1}_{m}^{\\!\mathrm{T}}\\!+\\!\mathbf{1}_{m}(\mathrm{vec}(X))^{\\!\mathrm{T}}\hat{L}^{\\!\mathrm{T}}),$

where $\mathcal{S}(\\!X\\!)\\!\in\\!K_{m}$ and
$\hat{L}\mathrm{vec}(\\!X\\!)\mathbf{1}_{m}^{\\!\mathrm{T}}\\!+\\!\mathbf{1}_{m}(\\!\mathrm{vec}(\\!X\\!)\\!)^{\\!\mathrm{T}}\\!\hat{L}^{\\!\mathrm{T}}\\!\in
S_{m}$. So the projection of $\mathcal{T}(X)$ onto $K_{m}$ is just
$\mathcal{S}(X)$. $\blacksquare$

Denote by $B_{1}$, $B_{2}$ and $B_{3}$ the base of $S_{m}^{0}$, $D_{m}$ and
$K_{m}$, respectively. By conclusion (i) of Proposition 3.1, we have
$\mathcal{T}(S_{m}^{0})\subset S_{m}^{0}$ and $\mathcal{T}(D_{m})\subset
S_{m}^{0}$. So

$\mathcal{T}[B_{1},\ B_{2},\ B_{3}]=[B_{1},\ B_{2},\
B_{3}]\left[\begin{array}[]{ccc}T_{11}&T_{12}&T_{13}\\\ &0&T_{23}\\\
&&T_{33}\end{array}\right].$ (11)

By the properties of the Lyapunov mapping $\mathcal{S}(\cdot)$ [2], [3], we
get

$\mathcal{S}[B_{1},\ B_{2},\ B_{3}]=[B_{1},\ B_{2},\
B_{3}]\left[\begin{array}[]{ccc}S_{11}&S_{12}&0\\\ S_{21}&S_{22}&0\\\
&&S_{33}\end{array}\right].$ (12)

From conclusion (ii) of Proposition 3.1, it follows that $T_{33}=S_{33}$.
Fortunately, the eigenvalues of $S_{33}$, i.e. the spectrum of the Lyapunov
mapping $\mathcal{S}(\cdot)$ restricted on $K_{m}$, have been revealed.

######  Lemma 3.1.

[3] Assume that the eigenvalues of $L$ be $\lambda_{1}$, $\lambda_{2}$,
$\lambda_{3}$, $\dots$, $\lambda_{m}$. Then the eigenvalues of the Lyapunov
mapping $\mathcal{S}(\cdot)$ restricted on $K_{m}$ are
$\\{\lambda_{i}+\lambda_{j}\ |\ 1\leq i<j\leq m\\}$.

Therefore, if all the eigenvalues of $\mathcal{T}(\cdot)$ on
$\mathbf{R}^{m\times m}$ are obtained, then by (11) and Lemma 3.1 we can get
the eigenvalues of $T_{11}$, i.e. the eigenvalues of $\mathcal{T}(\cdot)$
restricted on $S_{m}^{0}$. So, the rest content of this paper is just pursing
the eigenvalues of $\mathcal{T}(\cdot)$ on $\mathbf{R}^{m\times m}$. To this
end, we use the isomorphic mapping vec($X$): $\mathbf{R}^{m\times
m}\rightarrow\mathbf{R}^{m^{2}}$ to deal with $\mathcal{T}(\cdot)$.

######  Lemma 3.2.

[13] Let $A$, $B$ and $C$ be $m\times n$, $n\times s$ and $s\times t$
matrices, respectively. Then

$\mathrm{vec}(ABC)=(C^{\mathrm{T}}\otimes A)\mathrm{vec}(B).$

Applying Lemma 3.2 to $\mathcal{T}(X)$ described by (10) yields

$\mathrm{vec}(\mathcal{T}(X))\\!=\\!(L\otimes I_{m}+I_{m}\otimes
L-\mathbf{1}_{m}\otimes\hat{L}-\hat{L}\otimes\mathbf{1}_{m})\mathrm{vec}(X).$

So the eigenvalues are just those of the matrix

$T=L\otimes I_{m}+I_{m}\otimes
L-\mathbf{1}_{m}\otimes\hat{L}-\hat{L}\otimes\mathbf{1}_{m}.$ (13)

Before the main results, we first give a lemma as follows:

######  Lemma 3.3.

Let $A,B,C\in\mathbf{R}^{n\times n}$ satisfy $AB=BC$ and $B^{2}=0$. Then the
eigenvalue set of $A+B$ is the same as that of $A$.

Let rank$B=r$. When $r=0$ the assertion of the lemma is obviously right. When
$r>0$, there is an $n$-by-$r$ matrix $M$ and a $r$-by-$n$ matrix $N$ such that

$B=MN,\ \ \mathrm{rank}M=\mathrm{rank}N=r.$ (14)

Since $M$ has a full column rank, one can construct a nonsingular matrix
$T=[M,T_{1}]$. Let $T^{-1}=[P^{\mathrm{T}},Q^{\mathrm{T}}]^{\mathrm{T}}$. Then

$\left[\begin{array}[]{cc}I_{r}&0\\\
0&I_{n-r}\end{array}\right]=T^{-1}T=\left[\begin{array}[]{cc}PM&PT_{1}\\\
QM&QT_{1}\end{array}\right],$ (15)

which implies that

$QM=0,\ \ QB=QMN=0.$ (16)

Since $N$ has a full row rank and $BMN=B^{2}=0$, we have that $BM=0$. Thus

$T^{-1}BT=\left[\begin{array}[]{c}P\\\ Q\end{array}\right]B[M\ \
T_{1}]=\left[\begin{array}[]{cc}0&PBT_{1}\\\ 0&0\end{array}\right].$ (17)

From $AMN=AB=BC$ and (16), it follows that

$QAM=QBCN^{\mathrm{T}}(NN^{\mathrm{T}})^{-1}=0.$ (18)

Thus

$T^{-1}AT=\left[\begin{array}[]{c}P\\\ Q\end{array}\right]A[M\ \
T_{1}]=\left[\begin{array}[]{cc}PAM&PAT_{1}\\\ 0&QAT_{1}\end{array}\right].$
(19)

By (17) and (19), the proof is complete. $\blacksquare$

######  Proposition 3.2.

Assume that the Laplacian matrix $L$ satisfies rank$(L)=m-1$ and the
eigenvalues of $L$ are $\lambda_{1}=0$, $\lambda_{2}$, $\lambda_{3}$, $\dots$,
$\lambda_{m}$. Then the characteristic polynomial of $\mathcal{T}(\cdot)$ is

$s^{m}\prod_{i=2}^{m}(s-\lambda_{i})\prod_{j=2}^{m}(s-2\lambda_{j})\\!\\!\\!\prod_{2\leq
q<p\leq m}\\!\\!\\!(s-\lambda_{p}-\lambda_{q})^{2}.$ (20)

Since rank$(L)=m-1$ and $L\mathbf{1}_{m}=0$, there exists a nonsingular matrix
$P=[\mathbf{1}_{m},P_{1}]$ which results the Jordan canonical form of $L$ as
follows:

$J\\!=\\!P^{-1}\\!LP\\!=\\!\left[\\!\begin{array}[]{cc}0&\\\
&\tilde{J}\end{array}\\!\right]\\!\\!=\\!\\!\left[\\!\\!\begin{array}[]{c|cccc}0&\\\
\hline\cr&\lambda_{2}&b_{2}&&\\\ &&\ddots&\ddots&\\\ &&&\ \lambda_{m-1}&\ \
b_{m-1}\\\ &&&&\lambda_{m}\end{array}\\!\\!\right],$ (21)

where every $\lambda_{i}$ is nonzero and each $b_{i}$ is one or zero. It is
easy to check that

$\displaystyle(P^{-1}\otimes P^{-1})(L\otimes I_{m}+I_{m}\otimes
L-\mathbf{1}_{m}\otimes\hat{L})(P\otimes P)$ (28) $\displaystyle=$
$\displaystyle J\otimes I_{m}+I_{m}\otimes J-\delta_{1}\otimes
P^{-1}\hat{L}(P\otimes P)$ $\displaystyle=$ $\displaystyle J\otimes
I_{m}+I_{m}\otimes J-\delta_{1}\otimes P^{-1}\hat{L}[\mathbf{1}_{m}\otimes P,\
\ P_{1}\otimes P]$ $\displaystyle=$ $\displaystyle J\otimes I_{m}+I_{m}\otimes
J-\delta_{1}\otimes[P^{-1}LP,\ \ P^{-1}\hat{L}(P_{1}\otimes P)]$
$\displaystyle=$ $\displaystyle J\otimes I_{m}+I_{m}\otimes
J-\delta_{1}\otimes[J,\ \ P^{-1}\hat{L}(P_{1}\otimes P)]$ $\displaystyle=$
$\displaystyle\\!\\!\left[\\!\\!\begin{array}[]{cc}0_{m\\!\times\\!m}&\\\
&\\!\\!\\!\tilde{J}\otimes
I_{m}\end{array}\\!\\!\right]\\!\\!+\\!\\!\left[\\!\begin{array}[]{cc}J&\\\
&I_{m-1}\\!\otimes\\!J\end{array}\\!\right]\\!\\!-\\!\\!\left[\\!\\!\begin{array}[]{cc}J&\
\ P^{-1}\\!\hat{L}(P_{1}\\!\otimes\\!P)\\\ 0&0\end{array}\\!\\!\right]$
$\displaystyle=$ $\displaystyle\\!\left[\\!\begin{array}[]{cc}0_{m\times
m}&-P^{-1}\hat{L}(P_{1}\otimes P)\\\ &\ \ \tilde{J}\otimes
I_{m}+I_{m-1}\\!\otimes J\end{array}\\!\right].$ (31)

From (28), it follows that the characteristic polynomial of $L\otimes
I_{m}+I_{m}\otimes L-\mathbf{1}_{m}\otimes\hat{L}$ is just (20).

In the following, we use Lemma 3.3 to complete the proof. Let

$\displaystyle~{}~{}A$ $\displaystyle=$ $\displaystyle L\otimes
I_{m}+\\!I_{m}\otimes L-\mathbf{1}_{m}\otimes\hat{L},$ $\displaystyle~{}~{}B$
$\displaystyle=$
$\displaystyle\hat{L}\otimes\mathbf{1}_{m}=\left[\begin{array}[]{cccc}\mathbf{1}_{m}L_{1}&\\\
&\mathbf{1}_{m}L_{2}&\\\ &&\ddots\\\
&&&\mathbf{1}_{m}L_{m}\end{array}\right].$ (36)

From $L\mathbf{1}_{m}=0$, it follows that $B^{2}=0$, $(I_{m}\otimes L)B=0$ and
$(\mathbf{1}_{m}\otimes\hat{L})B=0.$ Now, we construct a matrix $C$ such that
$(L\otimes I_{m})B=BC$, which is equivalent to

$l_{ij}\mathbf{1}_{m}L_{j}=\mathbf{1}_{m}L_{i}C_{ij},\ \ \forall\ 1\leq
i,j\leq m.$ (37)

When $L_{i}=0$, we have $l_{ij}=0$, which means that $C_{ij}$ can be any
$m$-by-$m$ matrix. When $L_{i}\neq 0$, it is easy to check that
$C_{ij}=l_{ij}L_{i}^{\mathrm{T}}(L_{i}L_{i}^{\mathrm{T}})^{-1}L_{j}$ satisfies
(37). So, there exists a matrix $C$ such that $AB=BC$. Therefore, by Lemma
3.3, $A-B$ and $A$ have the same characteristic polynomial. The proof is
complete. $\blacksquare$ By Proposition 3.2, we have obtained all the
eigenvalues of the linear transformation $\mathcal{T}(\cdot)$. Now we can
determine the exact exponential synchronization rate.

######  Theorem 3.1.

Consider the high-dimensional Kuramoto model (3) limited on the unit sphere
$\mathbf{S}^{n-1}$. If the interconnecting digraph with the adjacency matrix
$(a_{ij})$ has a spanning tree, then the exponential synchronization rate is
exactly the smallest non-zero real part of the eigenvalues of Laplacian matrix
$L$.

Since the interconnecting digraph with the adjacency matrix $(a_{ij})$ has a
spanning tree, the Laplacian matrix $L$ satisfies rank$L=m-1$ [17]. So by
Proposition 3.2, we get the characteristic polynomial of $\mathcal{T}(\cdot)$
as shown in (20), where $\lambda_{1}=0$, $\lambda_{2}$, $\dots$, $\lambda_{m}$
are the eigenvalues of $L$. By Lemma 3.2 and $\lambda_{1}=0$, we get the the
characteristic polynomial of $T_{33}$ in (11) as follows:

$\prod_{i=2}^{m}(s-\lambda_{i})\\!\\!\prod_{2\leq q<p\leq
m}\\!\\!\\!(s-\lambda_{p}-\lambda_{q}).$ (38)

From (11), (20) and (38), it follows that the characteristic polynomial of
$T_{11}$ is

$\prod_{j=2}^{m}(s-2\lambda_{j})\\!\\!\prod_{2\leq q<p\leq
m}\\!\\!\\!(s-\lambda_{p}-\lambda_{q}),$ (39)

Without loss of generality, we assume that $\lambda_{2}$ has the smallest
nonzero real part of the eigenvalues of $L$. Therefore, by (39) we see that
the decay rate of $E(t)$ is $2\mathrm{Re}(\lambda_{2})$, which implies that
$|e_{ij}(t)|\leq c\mathrm{e}^{-2\mathrm{Re}(\lambda_{2})t}$ for a constant
$c>0$. Considering
$|e_{ij}(t)|=e_{ij}(t)=\frac{1}{2}||r_{i}(t)-r_{j}(t)||^{2}$, we conclude that
the exponential synchronization rate is $\mathrm{Re}(\lambda_{2})$.
$\blacksquare$

######  Remark 3.1.

The graph condition that the digraph has a directed spanning tree is the
weakest condition for synchronization. Actually, if the digraph has no a
directed spanning tree, there exist at least two independent strongly
connected components. Since there is no any information interaction between
two independent strongly connected components, it is impossible to achieve
synchronization for all the possible initial states.

##  4 Conclusions

For the high-dimensional Kuramoto model with identical oscillators and a
general digraph admitting a spanning tree, the exponential synchronization
rate has been accurately determined by using the matrix Riccati differential
equation of the synchronization error dynamics. In our future work, we will
try to generalize the main results and method to non-identical Kuramoto
oscillators and some other generalized high-dimensional Kuramoto models.

## References

  * [1] J. C. Bronski, T. E. Carty, and S. E. Simpson.  A matrix-valued Kuramoto model.  Journal of Statistical Physics, 178(2):595–624, 2020.
  * [2] D. Cheng.  On Lyapunov mapping and its applications.  Communications in Information and Systems, 1(3):255–272, 2001.
  * [3] D. Cheng and W. Xue.  Some properties of Lyapunov mapping.  IFAC Proceedings Volumes, 32(2):3724–3729, 1999.  14th IFAC World Congress 1999, Beijing, China, 5-9 July.
  * [4] D. Chi, S.-H. Choi, and S.-Y. Ha.  Emergent behaviors of a holonomic particle system on a sphere.  Journal of Mathematical Physics, 55(5):052703, 2014.
  * [5] S.-H. Choi and S.-Y. Ha.  Complete entrainment of Lohe oscillators under attractive and repulsive couplings.  SIAM Journal on Applied Dynamical Systems, 13(4):1417–1441, 2014\.
  * [6] N. Chopra and M. W. Spong.  On exponential synchronization of Kuramoto oscillators.  IEEE Transactions on Automatic Control, 54(2):353–357, 2009.
  * [7] J.-G. Dong and X. Xue.  Synchronization analysis of Kuramoto oscillators.  Communications in mathematical sciences, 11(2):465–480, 2013.
  * [8] Y. Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In H. Araki, editor, International Symposium on Mathematical Problems in Theoretical Physics, pages 420–422, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.
  * [9] C. Lageman and Z. Sun.  Consensus on spheres: Convergence analysis and perturbation theory.  In 2016 IEEE 55th Conference on Decision and Control (CDC), pages 19–24, 2016.
  * [10] W. Li and M. W. Spong.  Unified cooperative control of multiple agents on a sphere for different spherical patterns.  IEEE Transactions on Automatic Control, 59(5):1283–1289, 2014.
  * [11] M. A. Lohe.  Non-Abelian Kuramoto models and synchronization.  Journal of Physics A: Mathematical $\&$ Theoretical, 42(39):395101, 2009.
  * [12] M. A. Lohe.  Quantum synchronization over quantum networks.  Journal of Physics A: Mathematical $\&$ Theoretical, 43(46):465301, 2010.
  * [13] J. Magnus and H. Neudecker.  Matrix Differential Calculus With Applications in Statistics and Econometrics.  New York: John Wiley, 1988.
  * [14] J. Markdahl, J. Thunberg, and J. Goncalves.  Almost global consensus on the $n$-sphere.  IEEE Transactions on Automatic Control, 63(6):1664–1675, 2018.
  * [15] J. Markdahl, J. Thunberg, and J. Goncalves.  High-dimensional Kuramoto models on Stiefel manifolds synchronize complex networks almost globally.  Automatica, 113:108736, 2020.
  * [16] R. Olfati-Saber.  Swarms on sphere: A programmable swarm with synchronous behaviors like oscillator networks.  In Proceedings of the 45th IEEE Conference on Decision and Control, pages 5060–5066, 2006.
  * [17] W. Ren and R. W. Beard.  Consensus seeking in multiagent systems under dynamically changing interaction topologies.  IEEE Transactions on Automatic Control, 50(5):655–661, 2005.
  * [18] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths.  The Kuramoto model in complex networks.  Physics Reports, 610:1–98, 2016.
  * [19] A. Sarlette and R. Sepulchre.  Consensus optimization on manifolds.  SIAM Journal on Control and Optimization, 48(1):56–76, 2009.
  * [20] D. Tong, P. Rao, Q. Chen, M. J. Ogorzalek, and X. Li.  Exponential synchronization and phase locking of a multilayer Kuramoto-oscillator system with a pacemaker.  Neurocomputing, 308:129–137, 2018.
  * [21] Y. Wang and F. J. Doyle.  Exponential synchronization rate of Kuramoto oscillators in the presence of a pacemaker.  IEEE Transaction on Automatic Control, 58(4):989–994, 2013.
  * [22] J. Zhang and J. Zhu.  Exponential synchronization of the high-dimensional Kuramoto model with identical oscillators under digraphs.  Automatica, 102:122–128, 2019.
  * [23] J. Zhang, J. Zhu, and C. Qian.  On equilibria and consensus of the Lohe model with identical oscillators.  SIAM Journal on Applied Dynamical Systems, 17(2):1716–1741, 2018\.
  * [24] J. Zhu.  Synchronization of Kuramoto model in a high-dimensional linear space.  Physics Letters A, 377(41):2939–2943, 2013.
  * [25] J. Zhu.  High-dimensional Kuramoto model limited on smooth curved surfaces.  Physics Letters A, 378(18):1269–1280, 2014.