# Multivariate Wasserstein Functional Connectivity
for Autism Screening

###### Abstract

Most approaches to the estimation of brain functional connectivity from the
functional magnetic resonance imaging (fMRI) data rely on computing some
measure of statistical dependence, or more generally, a distance between
univariate representative time series of regions of interest (ROIs) consisting
of multiple voxels. However, summarizing a ROI’s multiple time series with its
mean or the first principal component (1PC) may result to the loss of
information as, for example, 1PC explains only a small fraction of variance of
the multivariate signal of the neuronal activity.

We propose to compare ROIs directly, without the use of representative time
series, defining a new measure of multivariate connectivity between ROIs, not
necessarily consisting of the same number of voxels, based on the Wasserstein
distance. We assess the proposed Wasserstein functional connectivity measure
on the autism screening task, demonstrating its superiority over commonly used
univariate and multivariate functional connectivity measures.

Index Terms—  functional magnetic resonance imaging, fMRI, functional
connectivity, optimal transport, clinical neuroscience

##  1 Introduction

Brain and its activity is a complex network consisting of neuronal
populations, coactivity of which have to be consistent across large
populations with individual differences and dependence on variables of
interest such as age, socioeconomic status, and presence or absence of various
neurodevelopmental and neurodegenerative deceases and conditions [1]. This
activity known as the functional connectivity (FC), opposed to the structural
connectivity which is the actual physical network of the neuronal populations
and their interconnections, is estimated from the resting state functional
magnetic resonance imaging (rs-fMRI).

Deceases and conditions are shown to be related to the increased, decreased,
abnormal or generally altered functional connectivity, with the autism
spectrum disorder (ASD) being prominent example. ASD is a condition related to
the brain development, which affects the individual’s social behavior. ASD
typically develops in infancy and early childhood and screened by the
behavioral tests which may be subjective. Therefore, it is important to bring
the quantitative, interpretable methods to screen autism, with fMRI being the
promising modality to address that purpose.

The common approach to estimate FC consist of dividing the raw fMRI image time
series into ROI groups defined by parcellation atlas. Next the representative
time series of each ROI are obtained by averaging or principal component
analysis (PCA), and some measure of dependence between each pair of
representatives, most commonly the Pearson correlation coefficient is
computed. We refer to this pipeline as the univariate approach, based on
univariate time series representatives, assuming the homogeneity of the time
series in each ROI.

Although as commonly raw time series are grouped into ROI based not on their
similarity, but according to atlases based on neuroanatomy or brain’s
cytoarchitecture, one could not expect time series homogeneity within a ROI or
single mode of their distribution. Recently introduced multivariate approach
[1, 2, 3, 4] addresses this issue, with the most popular test for multivariate
independence characterizing the relation strength is the distance correlation
[5].

In our work we propose and assess a new multivariate measure of the functional
connectivity based on the Wasserstein distance [6, 7] with the application to
ASD screening.

##  2 Materials and methods

###  2.1 Functional connectivity

The functional connectivity is estimated from fMRI data using some measure of
dependence between time series of each ROI, most often Pearson correlation
coefficient between the represetative time series of each ROI obtained by
averaging or taking the first principal component.

More formally, given $\\{X_{i}\\}_{i=1}^{N}$ random variables with equal
number of observations $t$, but of different dimensions
$\\{n_{i}\\}_{i=1}^{N}$ where $N$ is the number of ROIs, one seeks to estimate
some measure of statistial dependence $\rho_{ij}$ for each pair
$(X_{i},X_{j})$. If such estimating is done using $1$-dimensional summaries of
each $X_{i}$ we call that _univariate_ approach, otherwise if $X_{i}$ and
$X_{j}$ are compared directly we call that _multivariate_ approach.

Fig. 1: Left: Computation of the functional connectivity. In the more common
univariate approach a representative time series are computed for each ROI,
for example by averaging or PCA and are compared. In the multivariate approach
ROIs are compared directly, for example by distance correlation [5] or by the
Wasserstein distance as in our proposal. Right: Mean connectivity matrices
given by the Pearson and distance correlation and the Wasserstein distance.
For the visualization matrices were transformed to have ones on the diagonal
to facilitate the comparison.

####  2.1.1 Pearson correlation

Most often used similarity metric for comparing fMRI time series is _Pearson
correlation_ , which is defined for population as

$\rho(X,Y)=\frac{\mathrm{Cov}(X,Y)}{\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}}.$
(1)

Empirical Pearson correlation is given by

$\displaystyle r(\mathbf{x},\mathbf{y})$
$\displaystyle=\frac{\sum_{i=1}^{n}(x_{i}-\bar{\mathbf{x}})(y_{i}-\bar{\mathbf{y}})}{\sqrt{\sum_{i=1}^{n}x_{i}^{2}}\sqrt{\sum_{i=1}^{n}y_{i}^{2}}}$
(2)
$\displaystyle=\frac{\langle\mathbf{x},\mathbf{y}\rangle}{\|\mathbf{x}\|_{2}\|\mathbf{y}\|_{2}},$
(3)

with (3) valid for z-transformed vectors, standardized to the zero mean and
unit variance.

####  2.1.2 Distance correlation

Pearson’s correlation coefficient captures only linear association, with rank
correlation coefficients traditionally used to address this issue. The
relatively recently introduced distance correlation [5] in addition to being
able to capture nonlinear dependencies is a multivariate measure of
association between the groups of random variables of arbitrary, not of
necessarily equal size. The idea of distance correlation is based on measuring
the correlation of $t^{2}$ pairwise distances between realizations of random
variables $X\in\mathbb{R}^{n}$ and $Y\in\mathbb{R}^{m}$, with population
_distance correlation_ is defined analogously to the Pearson correlation

$dCor(X,Y)^{2}=\frac{\mathrm{dCov}^{2}(X,Y)}{\sqrt{\mathrm{dVar}^{2}(X)\mathrm{dVar}^{2}(Y)}}.$
(4)

For an observed random sample, given by the ROI matrices
$\mathbf{X},\mathbf{Y}=\\{(X_{k},Y_{k})\mid 1,\dots,n\\}$ from the joint
distribution of random vectors $X\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}^{q}$
define

$\displaystyle a_{kl}=|X_{k}-X_{l}|_{p},$
$\displaystyle\bar{a}_{k\cdot}=\frac{1}{n}\sum_{l=1}^{n}a_{kl},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ $\displaystyle\bar{a}_{\cdot
k}=\frac{1}{n}\sum_{k=1}^{n}a_{kl},$
$\displaystyle\bar{a}_{\cdot\cdot}=\frac{1}{n^{2}}\sum_{k,l=1}^{n}a_{kl}$

$k,l=1,\dots,n.$ to construct double centered pairwise distance matrix
$\mathbf{A}$ (and simirarly for $\mathbf{B}$) with entries given by

$(\mathbf{A})_{kl}=a_{kl}-\bar{a}_{k\cdot}-\bar{a}_{\cdot
k}+\bar{a}_{\cdot\cdot}.$ (5)

The empirical distance covariance $dCov_{n}(\mathbf{X},\mathbf{Y})$ is defined
by

$dCov_{n}^{2}(\mathbf{X},\mathbf{Y})=\frac{1}{n^{2}}\sum_{k,l=1}^{n}A_{kl}B_{kl}.$
(6)

The empirical distance variance $dVar_{n}(\mathbf{X})$ is defined by

$dVar_{n}^{2}(\mathbf{X})=\frac{1}{n}\sum_{k,l=1}^{n}A_{kl}^{2}.$ (7)

Finally, the empirical _distance correlation_
$dCor_{n}(\mathbf{X},\mathbf{Y})$ is defined

$dCor_{n}(\mathbf{X},\mathbf{Y})=\frac{dCov_{n}^{2}(\mathbf{X})}{\sqrt{dVar_{n}^{2}(\mathbf{X})dVar_{n}^{2}(\mathbf{Y})}},$
(8)

for $dVar_{n}^{2}(\mathbf{X})dVar_{n}^{2}(\mathbf{Y})>0$ and zero otherwise.

####  2.1.3 Wasserstein distance

Optimal transport considers comparing the measures over the domain $X$. Given
a ground metric $d:X\times X\rightarrow\mathbb{R}$ optimal transport equips
the space of measures $\mathcal{P}(X)$ with a metric referred to as the
_Wasserstein distance_ [6, 7], which for any $\mu,\nu\in\mathcal{P}(X)$ and
$p\geq 1$ is defined as [8]

$W_{p}^{p}(\mu,\nu)=\inf_{\pi\in\Pi(\mu,\nu)}\int_{X\times
X}d^{p}(x,y)d\pi(x,y)$ (9)

where $W_{p}^{p}$ denotes the $p$-th power of $W_{p}$ and $\Pi(\mu,\nu)$ is
the set of probability measures on the product space $X\times X$ whose
marginals coincide with $\mu$ and $\nu$; namely

$\small\Pi(\mu,\nu)=\\{\pi\in\mathcal{P}(X\times X)\mid
P_{1}\\#\pi=\mu,P_{2}\\#\pi=\nu\\}.$

Given ROI matrices $\mathbf{X}\in\mathbb{R}^{n\times t}$ and
$\mathbf{Y}\in\mathbb{R}^{m\times t}$, consider the discrete measures
$\mu_{\mathbf{X}},\nu_{\mathbf{Y}}\in\mathcal{P}(X)$ that can be written as
linear combinations
$\mu_{\mathbf{X}}=\sum_{i=1}^{n}a_{i}\delta_{\mathbf{x}_{i}}$ and
$\nu_{\mathbf{Y}}=\sum_{j=1}^{m}b_{j}\delta_{\mathbf{y}_{j}}$ of Dirac’s
deltas centred at a finite number $n$ and $m$ of points
$(\mathbf{x}_{i})_{i=1}^{n}\in\mathbb{R}^{t}$ and
$(\mathbf{y}_{j})_{j=1}^{m}\in\mathbb{R}^{t}$. In order for $\mu$ and $\nu$ to
be probabilities, set the vector weights
$\mathbf{a}=(1/n,\dots,1/n)^{T}\in\mathbb{R}^{n}$ and
$\mathbf{b}=(1/m,\dots,1/m)^{T}\in\mathbb{R}^{m}$.

We define the _Wasserstein functional connectivity (WFC)_ as the Wasserstein
distance between the two discrete measures $\mu_{\mathbf{X}}$ and
$\nu_{\mathbf{Y}}$ with corresponding weight vectors $\mathbf{a}$ and
$\mathbf{b}$ corresponds to

$W_{p}^{p}(\mu_{\mathbf{X}},\nu_{\mathbf{Y}})=\min_{\mathbf{T}\in\Pi(\mathbf{a},\mathbf{b})}\langle\mathbf{T},\mathbf{M}\rangle_{F},$
(10)

where $\mathbf{M}\in\mathbb{R}^{n\times m}$ is the cost matrix with entries
$(\mathbf{M})_{ij}=d^{p}(\mathbf{x}_{i},\mathbf{y}_{j})$ are set to the
$\ell^{p}$-distance $d^{p}(\cdot,\cdot)$,
$\langle\mathbf{T},\mathbf{M}\rangle_{F}$ is the elementwise Frobenius inner
product, and $\Pi(\mathbf{a},\mathbf{b})$ denotes the transportation polytope

$\Pi(\mathbf{a},\mathbf{b})=\\{\mathbf{T}\in\mathbb{R}^{n\times
m}_{+}:\mathbf{T}\mathbf{1}_{m}=\mathbf{a},\mathbf{T}^{T}\mathbf{1}_{n}=\mathbf{b}\\}.$

###  2.2 Dataset

Autism Brain Imaging Data Exchange (ABIDE) [9] is a collaboration of $16$
international imaging sites that have aggregated and are openly sharing
neuroimaging data. The dataset is composed of structural and resting state
functional MRI (rs-fMRI) data along with an extensive array of phenotypic
information. We used the preprocessed version111http://preprocessed-
connectomes-project.org/abide/ with $879$ subjects total, with $474$
individuals suffering from ASD and $405$ typical controls. Functional
preprocessing was performed using the Configurable Pipeline for the Analysis
of Connectomes (CPAC) including the band-pass filtering and global signal
regression.

###  2.3 Experiment

Functional connectivity matrices were computed using the Pearson correlation
between the ROI representative time series obtained by averaging, and the
distance correlation in the voxel space and the Wasserstein distance in the
time space for full ROIs. Being symmetric with the same values on the
diagonal, matrices were further vectorized by taking the upper half triangle
for each matrix.

Following the best practices in fMRI data analysis described in Dadi et al.
[10], we used logistic regression with L$1$ and L$2$ regularization with the
grid search for the regularization parameter on the $[0.2,5]$ range. The data
was split into the train and test subsets with the $80/20$ ratio, with the
random downsampling of the largest class performed to mitigate the class
imbalance. The classification accuracy mean and standard deviation in percents
for the best pair of the regularization norm and regularization parameter
value averaged over $20$ runs are reported in the Table 1.

Also we conducted the analysis of how considered connectivity measures are
correlated between themselves. For that, first, mean connectivity matrices
were computed for each measure, averaged for all subjects in the dataset.
Second, the absolute differences between mean connectivity matrices of the
ASD-affected and control patients were computed. Then, for each set of
matrices the correlation was computed for each pair in the set, summarizing
the values for each pair in the Table 2.

##  3 Results

| Univariate | Multivariate
---|---|---
| Pearson r | dCor | Wasserstein
All sites | 64.09$\pm$0.00 | 62.95$\pm$0.00 | 65.54$\pm$0.00
Caltech |  56.02$\pm$8.84 |  54.37$\pm$7.33 |  63.82$\pm$4.31
KKI |  55.40$\pm$11.87 |  65.30$\pm$8.87 |  55.65$\pm$10.36
Leuven |  59.05$\pm$5.58 |  68.60$\pm$5.37 |  59.42$\pm$9.02
MaxMun |  55.88$\pm$11.17 |  80.14$\pm$4.23 |  63.18$\pm$7.10
NYU |  62.96$\pm$3.35 |  58.51$\pm$3.02 |  64.32$\pm$3.68
OHSU |  53.30$\pm$8.46 |  50.35$\pm$10.59 |  53.75$\pm$11.63
Olin |  68.60$\pm$7.74 |  60.35$\pm$13.45 |  63.60$\pm$9.52
Pitt |  70.67$\pm$3.70 |  66.92$\pm$5.89 |  73.04$\pm$5.41
SBL |  60.55$\pm$6.73 |  48.55$\pm$8.16 |  68.00$\pm$12.08
SDSU |  59.95$\pm$12.13 |  55.10$\pm$9.02 |  57.60$\pm$14.46
Stanford |  58.88$\pm$6.12 |  93.57$\pm$2.95 |  97.07$\pm$0.17
Trinity |  64.29$\pm$8.14 |  51.85$\pm$6.01 |  56.57$\pm$6.42
UCLA |  64.74$\pm$3.36 |  50.04$\pm$5.29 |  65.34$\pm$3.05
UM |  68.65$\pm$5.20 |  70.41$\pm$4.57 |  72.98$\pm$5.47
USM |  65.66$\pm$4.27 |  64.43$\pm$6.46 |  69.81$\pm$5.21
Yale |  69.43$\pm$9.02 |  56.35$\pm$6.26 |  64.76$\pm$8.04
Best | 4/16 | 3/16 | 9/16
Avg. rank | 2.06 | 2.50 | 1.44
Table 1: Comparison of the multivariate Wasserstein functional connectivity
measure to the univariate and multivariate measures, represented by the
Pearson and distance correlations respectively. Best classification accuracy,
% for the ASD/control classification task is reported for the $5$-fold cross-
validation with $20$ repeats on the ABIDE1 dataset. The three largest sites
are emphasized by the bold font.

By the experimental evaluation on the ABIDE dataset – the benchmark for autism
diagnosis, we obtained the following results:

  * •

WFC outperforms the univarate and multivariate methods (Table 1) for the whole
dataset composed of the data from all $16$ sites.

  * •

WFC is competitive sitewise, quantified by the best average rank, including
being ranked first on the $3$ largest sites – NYU (171 subjects), UM (82), USM
(61). For the Stanford site using the WFC metric leads to mean accuracy
increase of $38.19\%$ compared to the Pearson correlation.

  * •

WFC is more similar to dCor than to Pearson r, as shown with the visual
inspection (Fig. 1, right) and high correlation between them obtained for the
averaged connectivity matrices as reflected in the Table 2, despite WFC is
being defined in the time space, not the voxel space contrary to dCor.

‘

|  r/dDcor | r/WFC | dCor/WFC
---|---|---|---
Mean |  0.265 | -0.584 | -0.893
Difference |  0.009 | 0.116 | 0.398
Table 2: Correlations between mean connectivity matrices computed using
different connectivity measures and differences between ASD and control mean
matrices.

##  4 Discussion

Albeit shown to be beneficial for the downstream tasks, the main limitation of
the Wasserstein distance connectivity measure is the increased computational
complexity compared to both Pearson and distance correlations. This limitation
can be addressed with several complementary approaches.

First is the use of dimensionality reduction – we observed that applying PCA
on the time domain, taking $25$% principal components with the subsequent
computation of the Wasserstein distance does affect the value of WFC measure
less than $1$% for the worst case over the population. Second, one can replace
exact computation of the Wasserstein distance with the computation-efficient
approximations, making use of the entropic regularization [8] or tree metrics
[11]. Third, while it is not affecting computational complexity, the
computation of the connectivity matrix entries is independent of each other,
so is amenable to the distributed and parallel computation.

The key suggestion that Basti et al. [4] make is that multivariate
connectivity analysis benefits from increased functional inhomogeneity of time
series constituting a ROI. The corresponding analysis of the dependence of the
inhomogeneity and the performance of the classifiers using Wasserstein
distance as connectivity measure relative to the univariate methods based on
average and principal component analysis-based representatives is left for the
future work.

#### Conclusions

We have proposed the new Wasserstein distance-based measure of the
multivariate brain functional connectivity and demonstrated that it is
preferable over the commonly used univariate and multivariate connectivity
measures given by the Pearson and distance correlations respectively for the
task of autism diagnosis.

##  5 Compliance with Ethical Standards

This is a numerical simulation study conducted retrospectively using human
subject data made available by Di Martino et al. [9] for which no ethical
approval was required.

##  6 Acknowledgements

No funding was received for conducting this study. The authors have no
relevant financial or non-financial interests to disclose.

## References

  * [1] Linda Geerligs and Richard Henson,  “Functional connectivity and structural covariance between regions of interest can be measured more accurately using multivariate distance correlation,”  NeuroImage, vol. 135, pp. 16–31, 2016.
  * [2] Alexander Walther, Hamed Nili, Naveed Ejaz, et al.,  “Reliability of dissimilarity measures for multi-voxel pattern analysis,”  NeuroImage, vol. 137, pp. 188–200, 2016.
  * [3] Kwangsun Yoo, Monica Rosenberg, Stephanie Noble, et al.,  “Multivariate approaches improve the reliability and validity of functional connectivity and prediction of individual behaviors,”  NeuroImage, vol. 197, pp. 212–223, 2019.
  * [4] Alessio Basti, Hamed Nili, Olaf Hauk, Laura Marzetti, and Richard Henson,  “Multi-dimensional connectivity: a conceptual and mathematical review,”  NeuroImage, p. 117179, 2020.
  * [5] Gábor Székely, Maria Rizzo, and Nail Bakirov,  “Measuring and testing dependence by correlation of distances,”  The Annals of Statistics, vol. 35, no. 6, pp. 2769–2794, 2007.
  * [6] Leonid Kantorovich and Gennady Rubinshtein,  “On a space of totally additive functions,”  Vestnik of the St. Petersburg University: Mathematics, vol. 13, no. 7, pp. 52–59, 1958.
  * [7] Leonid Vaserstein,  “Markov processes over denumerable products of spaces, describing large systems of automata,”  Problemy Peredachi Informatsii, vol. 5, no. 3, pp. 64–72, 1969\.
  * [8] Gabriel Peyré, Marco Cuturi, et al.,  “Computational optimal transport,”  Foundations and Trends® in Machine Learning, vol. 11, no. 5-6, pp. 355–607, 2019.
  * [9] Adriana Di Martino, Chao-Gan Yan, Qingyang Li, et al.,  “The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism,”  Molecular Psychiatry, vol. 19, no. 6, pp. 659–667, 2014.
  * [10] Kamalaker Dadi, Gaël Varoquaux, et al.,  “Benchmarking functional connectome-based predictive models for resting-state fmri,”  NeuroImage, vol. 192, pp. 115–134, 2019.
  * [11] Tam Le, Makoto Yamada, Kenji Fukumizu, and Marco Cuturi,  “Tree-sliced variants of wasserstein distances,”  Advances in Neural Information Processing Systems, vol. 32, pp. 12304–12315, 2019.
  * [12] Nathalie Tzourio-Mazoyer et al.,  “Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain,”  NeuroImage, vol. 15, no. 1, pp. 273–289, 2002.
  * [13] Rémi Flamary et al.,  “Pot: Python optimal transport,”  Journal of Machine Learning Research, vol. 22, no. 78, pp. 1–8, 2021.