# Robust Hierarchical Patterns for identifying MDD patients: A Multisite Study

Dushyant Sahoo ∗Corresponding author at: Center for Biomedical Image Computing
and Analytic, Perelman School of Medicine, University of Pennsylvania,
Philadelphia, USA
E-mail address<EMAIL_ADDRESS>Center for Biomedical Image Computing and
Analytic, Perelman School of Medicine, University of Pennsylvania,
Philadelphia, USA  Mathilde Antoniades  Center for Biomedical Image Computing
and Analytic, Perelman School of Medicine, University of Pennsylvania,
Philadelphia, USA  Cynthia H.Y. Fu  School of Psychology, University of East
London, London, UK  Centre for Affective Disorders, Institute of Psychiatry,
Psychology and Neuroscience, King’s College London, London, UK  Christos
Davatzikos  Center for Biomedical Image Computing and Analytic, Perelman
School of Medicine, University of Pennsylvania, Philadelphia, USA

###### Abstract

Many supervised machine learning frameworks have been proposed for disease
classification using functional magnetic resonance imaging (fMRI) data,
producing important biomarkers. More recently, data pooling has flourished,
making the result generalizable across a large population. But, this success
depends on the population diversity and variability introduced due to the
pooling of the data that is not a primary research interest. Here, we look at
hierarchical Sparse Connectivity Patterns (hSCPs) as biomarkers for major
depressive disorder (MDD). We propose a novel model based on hSCPs to predict
MDD patients from functional connectivity matrices extracted from resting-
state fMRI data. Our model consists of three coupled terms. The first term
decomposes connectivity matrices into hierarchical low-rank sparse components
corresponding to synchronous patterns across the human brain. These components
are then combined via patient-specific weights capturing heterogeneity in the
data. The second term is a classification loss that uses the patient-specific
weights to classify MDD patients from healthy ones. Both of these terms are
combined with the third term, a robustness loss function to improve the
reproducibility of hSCPs. This reduces the variability introduced due to site
and population diversity (age and sex) on the predictive accuracy and pattern
stability in a large dataset pooled from five different sites. Our results
show the impact of diversity on prediction performance. Our model can reduce
diversity and improve the predictive and generalizing capability of the
components. Finally, our results show that our proposed model can robustly
identify clinically relevant patterns characteristic of MDD with high
reproducibility.

Keywords: Hierarchical Sparse Model, Generative-Discriminative, Domain
Adaption, Major Depressive Disorder, Functional Connectivity, fMRI

##  1 Introduction

MDD is one of the most widespread psychiatric disorders characterized by
persistent sadness, depressed mood, low self-esteem, sleep disturbances,
emotional changes, and loss of interest in pleasurable activities, causing
disruptions to daily life Belmaker and Agam (2008). In addition, MDD causes
more than $800,000$ deaths each year globally and is also the leading cause of
disability Otte et al. (2016). Understanding the mechanism of MDD is crucial
for effective diagnosis, treatment and prevention, and understanding the
functioning of the human brain in a depressive state compared to a healthy
one. Considering the breadth of symptoms of this disease, it follows that
disruptions within and across multiple brains systems and networks must be at
play. Indeed, much work has been done to understand the functional brain
changes associated with MDD. However, much remains unknown about the
pathophysiology of the disease and the rates of relapse and recurrence remain
high Mueller et al. (1999), Kessler (2012).

Previous studies have shown that MDD is associated with disruptions in
regional functional connectivity and abnormal functional integration of
distributed brain regions Greicius et al. (2007), Liu et al. (2013), Wu et al.
(2011), Zhu et al. (2012). More recent approaches using seed-based
connectivity, independent component analyses, network homogeneity and graph
theory for functional connectivity analyses have revealed similar findings-
disruptions in functional networks and in between functional networks across
specific region pairs in MDD. The brain networks exhibiting abnormal
interactions in MDD include the Default Mode Network (awareness of internal
states), Dorsal Attention Network (external awareness), Fronto-Parietal
Network (top-down regulation of attention and emotion), Salient Network
(salient events) and Affective Network (emotion processing) Ye et al. (2015),
Kaiser et al. (2015), Yan et al. (2019), Mulders et al. (2015), Iwabuchi et
al. (2015), Brakowski et al. (2017).

Many efforts have been made to build functional connectivity-based predictive
models for identifying network-based biomarkers of depression Craddock et al.
(2009), Zeng et al. (2012), Bhaumik et al. (2017), Rosa et al. (2015), Zhao et
al. (2020a). Majority of studies are based on multivariate pattern analysis of
functional connectivity. Dictionary learning is one such approach that can
model low rank, sparse group-level interpretable networks having high
predictive power Batmanghelich et al. (2011), Eavani et al. (2013). Eavani et
al. (2015) proposed sparse Connectivity Patterns (SCPs), which decompose the
connectivity matrix of each subject as a set of shared sparse interpretable
patterns and subject-specific information. They extended the method by adding
a discriminative term to classify young adults vs. children Eavani et al.
(2014). More recently, D’Souza et al. (2020) used SCPs based generative-
discriminative model to predict clinically relevant networks characteristic of
Autism Spectrum Disorder.

Accurately classifying case-control groups has proven challenging due to
inter-patient variability Benkarim et al. (2021). In addition, biomarkers and
brain patterns learned for prediction from small and homogeneous datasets are
difficult to replicate and poorly generalize to new cohorts. In order to
improve generalizability and replicability of brain patterns and evaluate a
hypothesis in multiple sites/settings, there have been several open-access
neuroimaging data-sharing initiatives Alexander et al. (2017), Biswal et al.
(2010), Casey et al. (2018), Di Martino et al. (2017). In these initiatives,
data is pooled from multiple sites to capture demographically diverse
populations therefore building heterogeneous datasets that are more likely to
reflect the wider population.

Surprisingly, multi-site studies have shown lower classification performance
and poor generalization to data from new cohorts compared to single-site
studies Arbabshirani et al. (2017), Chen et al. (2016), Nielsen et al. (2013).
Analyzing data from these initiatives poses an inherent challenge due to
variability introduced from the diverging backgrounds of the subjects and from
site differences in MRI scanner hardware and software Kostro et al. (2014),
Yamashita et al. (2019), Shinohara et al. (2017), Abraham et al. (2017),
Jovicich et al. (2016), Noble et al. (2017). The non-biological variability
introduced due to pooling of the data can affect the biomarkers or common
features extracted from fMRI data Yu et al. (2018), these include functional
connectivity Shinohara et al. (2017) and sparse hierarchical factors Sahoo and
Davatzikos (2021b). The non-biological variability can lead to decreased
statistical power, spurious results and difficulty in identifying robust
biomarkers depending on the task. In addition, the correlation between site
effects and biological predictors can lead to an incorrect inference of non-
biological differences as biological. Thus, many neuroimaging studies need to
develop robust models that remove the non-biological variance and extract
biologically relevant information.

We use hierarchical sparse connectivity patterns (hSCPs) Sahoo et al. (2020)
to extract the relation between various functional networks in the human
brain. The method proposes a correlation matrix decomposition strategy,
capturing sparse hierarchical components agnostic to the parcellation scheme.
The sparse representation helps identify shared co-activation patterns in all
the subjects while subject-specific coefficients model the dataset’s
heterogeneity. See Appendix A for extended discussion on related work. The
hSCP can find relations between various functional networks across various
populations, which can help in discovering functional incoherence caused due
to MDD.

This paper aimed to identify the hSCP biomarkers of MDD while reducing the
effects of diversity (age, sex, site) common to large pooled datasets. Our
work builds on the hSCP model by using the discriminative nature of the
subject-specific weights extracted from hSCPs to classify MDD via logistic
regression in a large multi-site study. To tackle the variability introduced
due to pooling of the datasets, we use robust to site hSCP (rshSCP) Sahoo and
Davatzikos (2021b). The rshSCP method captures linear site effects and uses
adversarial learning to reduce site effects in the subject-specific
coefficients. We extend robust to site hSCP and introduce discriminative
rshSCP (dis-rshSCP) to extract homogeneous components discriminative of MDD by
reducing heterogeneity introduced due to covariates (age, sex and site), which
are known to affect neuroimaging analysis Alfaro-Almagro et al. (2021), Duncan
and Northoff (2013). Experiments on real datasets show that our approach can
improve the split-sample and leave one site predictability power of the
components while retaining the reproducibility of the components, thus
capturing informative heterogeneity. The classification performance on unseen
data indicates the generalizability of the model. Our results demonstrate that
MDD is associated with increased and decreased representation in patterns
associated with various functional networks. The results demonstrate our
framework’s potential in identifying patient-predictive biomarkers of a MDD.

Outline: We start by reviewing hierarchial Sparse Connectivity Patterns
(hSCPs). Then, in Section 3, we present our method, extracting interpretable
hSCPs which are discriminative of MDD and are robust to covariates (age, sex
and site). In Section 4, we demonstrate that our method could extract hSCPs
with high reproducibility and prediction power. This is followed by a
discussion on the interpretability of the extracted patterns, limitations and
future work.

##  2 Preliminaries

We first introduce preliminaries that serve as a basis for our discussions.

#### Notations:

We follow the same notation as in Sahoo et al. (2020), Sahoo and Davatzikos
(2021a, b). $\mathbb{S}^{P\times P}_{++}$ denotes symmetric positive definite
matrices of size $P\times P$. $\mathbf{A}\geq 0$ denotes that all the elements
of matrix $\mathbf{A}$ are greater than or equal to $0$. $P\times P$ matrix
with all elements equal to one is denoted by $\mathbf{J}_{P}$. Identity matrix
of size $P\times P$ is denoted by $\mathbf{I}_{P}$ and
$\mathbf{A}\circ\mathbf{B}$ denotes element-wise product between two matrices
$\mathbf{A}$ and $\mathbf{B}$.

#### Problem Setup:

The fMRI data of the $i^{th}$ subject having $P$ regions and $T$ time points
is denoted by $\mathbf{X}^{i}\in\mathbb{R}^{P\times T}$ with total $N$ number
of subjects or participants. Let $\mathbf{\Theta}^{i}\in\mathbb{S}^{P\times
P}_{++}$ be the correlation matrix where $\mathbf{\Theta}^{i}_{m,o}$ stores
the correlation between time series of $m^{th}$ and $o^{th}$ node. Let there
be total $S$ sites in the multi-site data and $\mathcal{I}_{s}$ be the set
storing subjects from site $s$. Let $\mathbf{y}_{site}\in\mathbb{R}^{N\times
S}$ be site labels encoded in one-hot manner, and
$\mathbf{y}_{age}\in\mathbb{R}^{N}$, $\mathbf{y}_{sex}\in\mathbb{R}^{N}$ and
$\mathbf{y}_{mdd}\in\mathbb{R}^{N}$ be the vectors storing information about
age, sex and MDD label. We aim to extract set of hierarchical patterns
representative of depression using fMRI data with reduced variability
introduced due to age, sex and site.

$\cdots$ $\cdots$ $\mathbf{\Theta}^{n}$ $=$ $\mathbf{W}$ $*$
$\mathbf{\Lambda}^{n}$ $*$ $\mathbf{W}^{\top}$ $\top$ $\mathcal{D}$
$\mathbf{b}$ $\mathbf{y}_{mdd}$ $*$ $\xrightarrow{f(\mathcal{D},\mathbf{b})}$
Robustly Learning Components Classification Model Figure 1: A joint two level
modeling for connectivity analysis and prediction. First part is functional
data representation depicted in blue box. Here the correlation matrices are
decomposed into shared components stored in a basis matrix and subject
specific information. We learn this decomposition in a robust manner to reduce
variability due to site and demographics. To visualize each component, a
column of basis matrix is projected onto the brain. Second part is the
prediction of MDD patients depicted in green box.

###  2.1 Introduction to hierarchical Sparse Connectivity Patterns

A graphical summary of our model is presented in Fig. 1. The two inputs to our
model are the rs-fMRI correlation matrices (upper left) and the binary scores
depicting if a person is healthy or has MDD (lower right). The correlation
matrices are constructed from the time series data describing the similarity
using Pearson’s Correlation Coefficient between various nodes of the human
brain. The blue box in Fig 1 indicates the generative model estimating the
components. Here, we decompose correlation matrices into a set of components
capturing co-activation patterns common across the entire cohort and subject-
specific information capturing heterogeneity in the data representing the
strength of each component in each individual. We capture this information
while reducing the effects of demographics and site to improve the
generalizability of the components and predictability of the subject-specific
information. The green box indicates the discriminative model guiding the
components to represent MDD. Here, we use the information from the subject-
specific coefficients to predict MDD via a classification model for each
individual. $f()$ takes subject-specific information and a set of weights as
input and maps it to a binary value representing $1$ and $0$ for MDD and
healthy subjects.

Hierarchical Sparse Connectivity Patterns (hSCP) Sahoo et al. (2020) is an
extension of Sparse Connectivity Patterns (SCPs) Eavani et al. (2015) to find
hierarchical structures in the fMRI data to have multi-scale representation of
connectivity patterns in the human brain. It was first defined by Eavani et
al. (2015) to estimate sparse co-activated patterns in the human brain. As
seen in Fig. 1, we model $\mathbf{\Theta}^{n}$ using a shared representation
and a subject-specific term storing information about each component’s
strength. The $\mathbf{W}_{1}\in\mathbb{R}^{P\times k_{1}}$ stores $k_{1}$ is
a concatenation of $K$ components $\mathbf{w}_{l}\in\mathbb{R}^{P\times 1}$ ,
i.e.
$\mathbf{W}\coloneqq\left[\mathbf{w}_{1}\;\mathbf{w}_{2}\cdots\mathbf{w}_{k}\right]$,
where $k\leq P$. These components capture steady-state co-activation patterns
across regions in the brain shared across all the population. While components
are shared across all the subjects, their representation in each subject is
unique, captured by $\mathbf{\Lambda}_{1}^{n}$. $\mathbf{\Lambda}_{1}^{n}$ is
kept to be non-negative to ensure the positive semi-definite structure of the
correlation matrices $\mathbf{\Theta}^{n}$. The SCP representation is written
as:

$\displaystyle\mathbf{\Theta}^{n}\approx\mathbf{W}_{1}\mathbf{\Lambda}_{1}^{n}\mathbf{W}_{1}^{\top}.$

This formulation helps reduce the high dimensionality of the data by
projecting it into lower-dimensional space and providing subject-level
information. SCP can be extended to hSCP following the below equations:

$\displaystyle\mathbf{\Theta}^{n}$
$\displaystyle\approx\mathbf{W}_{1}\mathbf{\Lambda}_{1}^{n}\mathbf{W}_{1}^{\top},$
(1) $\displaystyle\mathbf{\Theta}^{n}$
$\displaystyle\approx\mathbf{W}_{1}\mathbf{W}_{2}\mathbf{\Lambda}_{1}^{n}\mathbf{W}_{2}^{\top}\mathbf{W}_{1}^{\top},$
$\displaystyle\vdots$ $\displaystyle\mathbf{\Theta}^{n}$
$\displaystyle\approx\mathbf{W}_{1}\mathbf{W}_{2}\ldots\mathbf{W}_{K}\mathbf{\Lambda}^{n}_{K}\mathbf{W}_{K}^{\top}\mathbf{W}_{K-1}^{\top}\ldots\mathbf{W}_{1}^{\top},$

where $K$ is the level of hierarchy in the formulation,
$\mathbf{\Lambda}^{n}_{k}$ is a non-negative diagonal matrix storing relative
contribution of each component for $n$th subject at $k$th level and
$P>k_{1}>\ldots>k_{K}$ to make the formulation low dimensional and forcing
each successive level in the hierarchy to have less number of components than
the previous level. In the hierarchical structure, at the bottom level we have
$\mathbf{W}_{1}\in\mathbb{R}^{P\times k_{1}}$ storing $k_{1}$ components.
These components are then linearly transformed at each successive level to get
lower-dimensional components than the previous level. $\mathbf{W}_{2}$,
$\mathbf{W}_{3}$, $\ldots$, $\mathbf{W}_{K}$ are used to transform $k_{1}$
dimensional components to $k_{2},k_{3},\ldots,k_{K}$ dimensional components
where $P>k_{1}>\ldots>k_{K}$. For example, let $\mathbf{W}_{1}$ be the
components at lower most level, then to extract components at level above
this, we multiply $\mathbf{W}_{1}$ by $\mathbf{W}_{2}$ to obtain $k_{2}$
dimensional components which are extracted using $k_{1}$ dimensional
components. All the elements of $\mathbf{W}_{2}$, $\mathbf{W}_{3}$, $\ldots$,
$\mathbf{W}_{K}$ are constrained to be greater than $0$ for interpretability.
We now define sets storing the correlation matrices, components and weights at
all level for notational brevity. Set of correlation matrices is denoted by
$\mathcal{C}={\\{\mathbf{\Theta}^{n}\mid n=1,\ldots,N\\}}$, set storing
information about shared components is denoted by
$\mathcal{W}={\\{\mathbf{W}_{r}\mid r=1,\ldots,K\\}}$ and the set storing
subject-specific information is denoted by
$\mathcal{D}={\\{\mathbf{\Lambda}_{r}^{n}\mid r=1,\ldots,K;n=1\ldots,N\\}}$
with $\mathbf{\Lambda}^{n}_{r}\geq 0$. The hierarchical components and
subject-specific information are estimated by performing alternating
minimization on the below optimization problem:

$\displaystyle\underset{\mathcal{W},\mathcal{D}}{\text{min}}$
$\displaystyle\sum_{n=1}^{N}\sum_{r=1}^{K}\|\mathbf{\Theta}^{n}-(\prod_{j=1}^{r}\mathbf{W}_{j})\mathbf{\Lambda}_{r}^{n}(\prod_{j=1}^{r}\mathbf{W}_{n})^{\top}\|_{F}^{2}$
(2) s.t.
$\displaystyle\|\mathbf{w}^{r}_{l}\|_{1}<\lambda_{r},\;l=1,\ldots,k_{r};\;r=1,\ldots,K$
$\displaystyle\|\mathbf{w}^{r}_{l}\|_{\infty}\leq
1,\;l=1,\ldots,k_{r};\;r=1,\ldots,K$ $\displaystyle\mathop{\rm
trace}\nolimits(\mathbf{\Lambda}_{r}^{n})=1,\;n=1,\ldots,N;\;r=1,\ldots,K$
$\displaystyle\mathbf{\Lambda}_{r}^{n}\geq 0,\;\;n=1,\ldots,N;\;r=1,\ldots,K$
$\displaystyle\mathbf{W}_{j}\geq 0,\;j=2,\ldots,K.$

The first constraint in the above optimization problem is $L_{1}$ which helps
in capturing compact and clinically informative representation by encouraging
sparsity in the components. $L_{1}$ regularizer will help in selecting a small
number of nonzero entries in $\mathbf{W}$ that explain the data.
Identifiability and reproducibility issues are solved by combining $L_{1}$,
$L_{\infty}$ and $\mathop{\rm trace}\nolimits$ constraints.
$\mathbf{\Lambda}_{r}^{n}\geq 0$ constraint ensures that the subject-specific
information is kept positive, and therefore the reconstructed matrix is
positive semi-definite. Again for notational brevity, we will use the set
notation for the above constraints by denoting
$\Omega_{\mathcal{W}}=\\{\mathbf{W}\mid\|\mathbf{w}^{r}_{l}\|_{1}<\lambda_{r},\;\|\mathbf{w}^{r}_{l}\|_{\infty}\leq
1,\mathbf{W}_{j}\geq 0,\;j=2,\ldots,K\\}$ and
$\Psi=\\{\mathbf{\Lambda}\mid\mathop{\rm
trace}\nolimits(\mathbf{\Lambda}_{r}^{n})=1,\mathbf{\Lambda}_{r}^{n}\geq
0\\}$. Refer to Sahoo et al. (2020) for a detailed description of the method.

###  2.2 Reducing site effects

In this section, we discuss framework proposed by Sahoo and Davatzikos (2021a)
to reduce site effects in the hSCP model. Their idea is to jointly model
sparse patterns and site information, and use adversarial learning to reduce
site effects. We first look at the case with one level in the hierarchy for
easier understanding. This will extended to multiple levels in the upcoming
section. Authors of the paper hypothesize that site and scanner information
for the training data can be represented by $\mathbf{V}\in\mathbb{R}^{P\times
P}$ space. For each site $s=1,\ldots,S$, let there be diagonal matrices
$\mathbf{U}^{s}\in\mathbb{R}^{P\times P}$ storing additive and
$\mathbf{Z}^{s}\in\mathbb{R}^{k\times k}$ storing scaled site-specific
information. Based on the above hypothesis, $\mathbf{\Theta}^{n}$ is
decomposed to jointly estimate the hSCPs, site and scanner information for
$n\in\mathcal{I}_{s}$ as:

$\displaystyle\mathbf{\Theta}^{n}\approx\mathbf{W}\mathbf{Z}^{s}\mathbf{\Lambda}^{n}\mathbf{W}^{\top}+\mathbf{U}^{s}\mathbf{V}.$
(3)

The method imposes $L_{1}$ sparsity constraint on $\mathbf{V}$ to prevent
overfitting.

##  3 Method

###  3.1 Robust to covariates hSCP

In the previous section, we saw how site and scanner effects could be stored
using $\mathbf{U}$ and $\mathbf{V}$. Sahoo and Davatzikos (2021a) added an
adversarial layer on top of it to reduce the predictive power of $\mathcal{D}$
storing subject-specific coefficients for predicting site. In addition to
reducing site effects in this paper, we also minimize variability due to age
and sex. For this, we train a model $F(\zeta,\mathcal{D})$ parameterized by
$\zeta$ with input $\mathcal{D}$ that return age
$\hat{\mathbf{y}}_{age}\in\mathbb{R}^{N}$, sex
$\hat{\mathbf{y}}_{sex}\in\mathbb{R}^{N}$ and site
$\hat{\mathbf{y}}_{site}\in\mathbb{R}^{N\times S}$ predictions. We train the
model by minimizing the below loss function
$\mathcal{L}_{1}(\zeta,\mathcal{D},\mathbf{y}_{site},\mathbf{y}_{age},\mathbf{y}_{sex})$:

$\displaystyle\mathcal{L}_{1}=\underbrace{\|\hat{\mathbf{y}}_{age}-{\mathbf{y}}_{age}\|_{2}^{2}}_{\begin{subarray}{c}\text{preserve
age}\\\
\text{information}\end{subarray}}-\sum_{s=1}^{S}\sum_{n=1}^{N}\underbrace{y^{n,s}_{site}\log\hat{y}^{n,s}_{site}}_{\begin{subarray}{c}\text{preserve
site}\\\
\text{information}\end{subarray}}-\sum_{n=1}^{N}\underbrace{\left(y^{n}_{sex}\log\hat{y}^{n}_{sex}+(1-y^{n}_{sex})\log(1-\hat{y}^{n}_{sex})\right)}_{\begin{subarray}{c}\text{preserve
sex}\\\ \text{information}\end{subarray}}.$ (4)

Let $\zeta^{*}=\underset{\zeta}{\text{arg min}}\;\mathcal{L}_{1}$ be optimum
value which minimizes $\mathcal{L}_{1}$. The above problem is the Multi-task
learning (MTL) Caruana (1997), Zhang and Yang (2017) problem to learn multiple
correlated tasks at the same time. This formulation helps improve the
performance of each and reduces the need to introduce multiple models to solve
individual tasks. We use the direct sum approach to combine different
objectives by directly minimizing the sum of all losses of different tasks, a
common practice in multi-task learning. A hard parameter sharing strategy is
used to model network architecture. In this strategy, parameters are shared by
the bottom layers among all tasks, while top layers are selected to be task-
specific, helping with robustness against overfitting Ruder (2017). It is a
commonly used method for designing deep learning models in the literature Long
et al. (2015), Ruder et al. (2019), Sener and Koltun (2018). Details about the
architecture are given in Appendix C.

Using the loss function $\mathcal{L}_{1}$, we modify $\mathbf{\Lambda}^{n}$
such that its predictability power to predict site, age and sex reduces. This
can be achieved by maximizing $\mathcal{L}_{1}$ loss with respect to
$\mathbf{\Lambda}^{n}$. Note that we are trying to solve two problems with one
loss function, first is finding optimal $\zeta$ which minimizes
$\mathcal{L}_{1}$ and second is finding optimal $\mathbf{\Lambda}^{n}$ which
maximizes $\mathcal{L}_{1}$. This will result in a minimax game, where the
$\zeta$ is learned to minimize the cross-entropy and regression loss, and
$\mathbf{\Lambda}^{n}$ is adjusted to maximize the loss. The minimax
optimization problem can be written as:

$\displaystyle\underset{\zeta}{\text{max}}$
$\displaystyle\underset{\mathbf{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z}}{\text{min}}\quad$
$\displaystyle\sum_{s=1}^{S}\sum_{n\in\mathcal{I}_{s}}\|\mathbf{\Theta}^{n}-\mathbf{W}\mathbf{Z}^{s}\mathbf{\Lambda}^{n}\mathbf{W}^{\top}-\mathbf{U}^{s}\mathbf{V}\|_{F}^{2}-\gamma_{1}\mathcal{L}_{1}(\zeta,\mathcal{D},\mathbf{y}_{site},\mathbf{y}_{age},\mathbf{y}_{sex})$
(5) $\displaystyle\hskip 10.00002pts.t.$
$\displaystyle\mathbf{W}\in\Omega,\quad\mathcal{D}\in\Psi,\quad\|\mathbf{v}_{p}\|_{1}<\mu,\;p=1,\ldots,P,$

where $\mathcal{U}=\\{\mathbf{U}_{s}|s=1,\ldots,S\\}$,
$\mathcal{Z}=\\{\mathbf{Z}_{s}|s=1,\ldots,S\\}$ and $\mathbf{v}_{p}$ is the
$p$th column of $\mathbf{V}$. Here
$\|\mathbf{\Theta}^{n}-\mathbf{W}\mathbf{Z}^{s}\mathbf{\Lambda}^{n}\mathbf{W}^{\top}-\mathbf{U}^{s}\mathbf{V}\|_{F}^{2}$
is the total error in the representation of the $S$ subjects and
$\mathcal{L}_{1}$ is the robustness loss, and $\gamma_{1}$ is the tradeoff
between representation learning and robustness.

###  3.2 Joint modeling of MDD scores

The aim of the paper is to learn components which are representative of MDD.
For this, we build discriminative rshSCP (dis-rshSCP) to use subject-specific
information $\mathbf{\Lambda}^{n}_{r}$ to predict whether subject $n$ has MDD
or not. Use this model we will subject-specific information which is most
predictive of MDD and will give components corresponding to that. We model MDD
information using logistic regression framework with parameter
$\mathbf{b}\in\mathbb{R}^{f}$, where $f=\sum_{r=1}^{K}k_{r}$, subject specific
information of all the levels combined
$\mathbf{t}^{n}=diag[\mathbf{\Lambda_{1}^{n}},\ldots,\mathbf{\Lambda_{K}^{n}}]$
and the loss function defined below:

$\displaystyle\mathcal{L}_{2}(\mathbf{b},\mathcal{D},\mathbf{y}_{mdd})=-\sum_{n=1}^{N}\left[y^{n}_{mdd}\log\left(\frac{1}{1+\exp{(-\mathbf{b}^{\top}\mathbf{t}^{n})}}\right)+(1-y^{n}_{mdd})\log\left(\frac{\exp{(-\mathbf{b}^{\top}\mathbf{t}^{n})}}{1+\exp{(-\mathbf{b}^{\top}\mathbf{t}^{n})}}\right)\right].$
(6)

We minimize cross entropy loss
$\mathcal{L}_{2}(\mathbf{b},\mathcal{D},\mathbf{y}_{mdd})$ with $\mathbf{b}$
as the parameters to be estimated. Now, let the generative loss be

$\displaystyle
G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})=\sum_{s=1}^{S}\sum_{n\in\mathcal{I}_{s}}\sum_{r=1}^{K}\|\mathbf{\Theta}^{n}-(\prod_{j=1}^{r}\mathbf{W}_{j})\mathbf{Z}^{s}_{r}\mathbf{\Lambda}_{r}^{n}(\prod_{j=1}^{r}\mathbf{W}_{n})^{\top}-\mathbf{U}^{s}_{r}\mathbf{V}_{r}\|_{F}^{2},$
(7)

which models the components and site information, then the joint optimization
problem can be written as:

$\displaystyle\underset{\zeta}{\text{max}}$
$\displaystyle\underset{\mathcal{W},\mathcal{D},\mathcal{U},\mathcal{Z},\mathbf{V},\mathbf{b}}{\text{min}}\quad$
$\displaystyle\underbrace{G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})}_{\begin{subarray}{c}\text{learn
subject and}\\\ \text{site
information}\end{subarray}}-\gamma_{1}\underbrace{\mathcal{L}_{1}(\zeta,\mathcal{D},\mathbf{y}_{site},\mathbf{y}_{age},\mathbf{y}_{sex})}_{\begin{subarray}{c}\text{reduce
age, sex and}\\\ \text{site
information}\end{subarray}}+\gamma_{2}\underbrace{\mathcal{L}_{2}(\mathbf{b},\mathcal{D},\mathbf{y}_{mdd})}_{\begin{subarray}{c}\text{preserve
MDD}\\\ \text{information}\end{subarray}}$ (8) $\displaystyle\hskip
10.00002pts.t.$
$\displaystyle\mathcal{W}\in\Omega,\quad\mathcal{D}\in\Psi,\quad\|\mathbf{v}_{p}\|_{1}<\mu,\;p=1,\ldots,P.$

Here, in addition to representation learning and robustness loss, we also have
prediction error $\mathcal{L}_{2}$. $\gamma_{1}$ and $\gamma_{2}$ are the
trade-offs between representation learning, robustness, and prediction.

###  3.3 Prediction on unseen data

To estimate $\mathbf{\Lambda}$ for a new subject, we first solve the
optimization problem in equation 5 to estimate $\mathbf{W}$ computed from the
training data. The estimation of the coefficients of unseen subjects are then
estimated by solving the below minimax problem where $\mathbf{W}$ is computed
from the training data:

$\displaystyle\underset{\zeta}{\text{max}}$
$\displaystyle\underset{\mathcal{D},\mathcal{U},\mathcal{Z},\mathbf{V}}{\text{min}}\quad$
$\displaystyle
G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})-\gamma_{1}\mathcal{L}_{1}(\zeta,\mathcal{D},\mathbf{y}_{site},\mathbf{y}_{age},\mathbf{y}_{sex})$
(9) $\displaystyle\hskip 10.00002pts.t.$
$\displaystyle\mathcal{D}\in\Psi,\quad\|\mathbf{v}_{p}\|_{1}<\mu,\;p=1,\ldots,P.$

The estimate for the MDD information for the test subject $n$ is given by:

$y_{mdd}^{n}=\begin{cases*}1&\text{if
$\frac{1}{1+\exp{(-\mathbf{b}^{\top}\mathbf{t}^{n})}}\geq 0.5$ }\\\
0&\text{otherwise},\end{cases*}$

where $\mathbf{b}$ and $\mathbf{t}$ is estimated from solving equation 8 and 9
respectively.

The optimization problems defined in 8 and above 9 are non-convex problems. We
use alternating minimization to solve, Appendix B gives the detail of the
complete algorithm and the optimization procedure. We see during convergence
that the reconstructions loss is a little higher if we hadn’t introduced
adversarial and discriminative losses. Here, we can expect a trade-off between
finding components with optimal reconstruction, adversarial and discriminative
loss that depends on $\gamma_{1}$ and $\gamma_{2}$. The optimal value of these
hyperparameters is selected using cross-validation, which we explain in the
Section 5.

##  4 Materials

###  4.1 Participants

Five worldwide study samples totaling $1657$ participants, including $733$
with MDD and $924$ healthy controls (HC) contributed T1-weighted structural
scans and resting-state fMRI data (rs-fMRI) to this study. The included
cohorts combine data from the following studies: EMBARC ($4$ centers across
the United States of America, Trivedi et al. (2016)), University of Oxford
(United Kingdom, Godlewska et al. (2014, 2018)), Sichuan University Cohort
(China, Zhao et al. (2020b)) and STRADL (United Kingdom, Navrady et al.
(2018), Stolicyn et al. (2020)). Patient and controls were on average $68.6\%$
($50-73.8\%$) and $54.5\%$ ($53-70\%$) female. The mean age across samples was
$47$ ($18-78$) years in patients and $61$ ($16-84$) years for controls. All
patients in EMBARC, Oxford and SCU were medication-free, and $15$ in SNAP and
$170$ in STRADL were medicated at the time of scanning and had a primary
diagnosis of MDD that was a first episode or recurrent. MDD diagnosis was
based on standardized diagnostic criteria: DSM-IV (Oxford) and DSM-IV-TR
(EMBARC, Stanford, STRADL, SCU) Frances et al. (1995), First et al. (2004).
Table 1 summarizes the number of healthy and MDD participants in each site
with their age and sex distribution.

Site | $N_{\textnormal{HC}}$ | $N_{\textnormal{MDD}}$
---|---|---
| Number | $\%\text{of F}$ | Age (y) | Number | $\%\text{of F}$ | Age (y)
EMBARC-CU | $12$ | $60.0$ | $[18,54](34)$ | $77$ | $68.3$ | $[18,64](30)$
EMBARC-MG | $10$ | $66.6$ | $[18,65](28)$ | $52$ | $56.2$ | $[18,64](28.5)$
EMBARC-TX | $11$ | $50.0$ | $[23,57](26)$ | $97$ | $67.7$ | $[19,65](44)$
EMBARC-UM | $10$ | $70.0$ | $[23,62](41.4)$ | $59$ | $71.1$ | $[18,65](31)$
Oxford | $31$ | $58.0$ | $[19,58](28)$ | $39$ | $60.5$ | $[20,61](27)$
SCU | $40$ | $55.0$ | $[16,57](26.5)$ | $30$ | $50.0$ | $[18,60](30.5)$
SNAP | $55$ | $65.3$ | $[19,58](28.8)$ | $55$ | $63.6$ | $[20,56](28.2)$
STRADL | $755$ | $53.1$ | $[26,84](62)$ | $324$ | $73.8$ | $[26,78](60)$
Total | $924$ | $54.5$ | $[16,84](61)$ | $733$ | $68.6$ | $[18,78](47)$
Table 1: Demographic characteristics of participants.

###  4.2 Data Preprocessing

We used FMRIB Software Jenkinson et al. (2012) and applied UK Biobank pipeline
as a part of initial pre-processing. The data pre-processing is the same as in
Sahoo and Davatzikos (2021b). The pre-processing included the removal of the
first five volumes, FSL’s MCFLIRT Jenkinson et al. (2012) for head movement
correction, global 4D mean intensity normalization, and temporal high-pass
filtering ($>0.01$ Hz). To remove structured artifacts, we used FMRIB’s ICA-
based Xnoiseifier FIX Salimi-Khorshidi et al. (2014), Griffanti et al. (2014).
FLIRT was then used to co-register functional images to the T1 image with
Boundary Based Registration (BBR) as the cost function and FSL’s FNIRT (non-
linear registration) was used to register T1-weighted images to the MNI152
template. We used Independent Components Analysis (ICA) maps Smith et al.
(2014) as our parcellation scheme to project the data into a $100$ dimensional
space i.e., each subject had $100$ nodes after applying the parcellation
scheme. Group Information Guided ICA (GIGICA) is used to extract one time
series per node by mapping the group independent component spatial maps onto
each subject’s time series Du and Fan (2013).

##  5 Experiment

###  5.1 Evaluating predictive performance

We use two different strategies to evaluate the performance of dis-rshSCP. In
the first strategy, we compare five-fold cross-validation accuracy. The goal
here is to check how well our model is able to estimate $\mathbf{\Lambda}$ to
classify MDD vs. healthy people. We train the model on $80\%$ training set for
each fold and test on the remaining. Here training is referred to solving
equation 8 to estimate
$\mathcal{W},\mathcal{D},\mathcal{U},\mathcal{Z},\mathbf{V}$ and $\mathbf{b}$.
During the test, we fix $\mathbf{b}$, estimate the rest of the parameters and
evaluate the performance of the estimated parameters using $\mathbf{b}$
estimated from training. We use stratification to ensure each fold is
representative of all strata (age, sex, site, and MDD) of the data. This is
performed to ensure each class is approximately equally represented across
each test and training fold. In the second strategy, we compare leave one site
out accuracy. We train the model on all sites except one and test the model on
the data from the remaining site. We compare the performance of 4 different
versions of the model: 1) vanilla hSCP without reduction of any covariates, 2)
dem-rshCP model with reduction of effects of age and sex, 3) site-rshSCP model
with a reduction in site effects and 4) com-rshSCP model with a reduction in
age, sex and site effects. We fix $k_{2}=4$ for better interpretability and
based on previous experiments of Sahoo et al. (2020) and Sahoo and Davatzikos
(2021b), and find optimal value of $k_{1}$ from the set $\\{5,10,15,20\\}$.
Note that even if we select a large value of $k_{2}$, only a few of those
components will be used to predict MDD and our experiments show that it is
less than $4$. Optimal value of hyperparameters $\mu$, $\gamma_{1}$,
$\gamma_{2}$ and $\lambda_{1}$ are selected from $[0.1,1,5]$, $[.5,1,5]$,
$[0.1,1,5]$ and $[0.1,1,10]$, respectively.

Method | $k_{1}=5$ | $k_{1}=10$ | $k_{1}=15$ | $k_{1}=20$
---|---|---|---|---
hSCP | $0.569\pm 0.024$ | $0.571\pm 0.022$ | $0.576\pm 0.022$ | $0.578\pm 0.019$
dem-rshSCP | $0.618\pm 0.021$ | $0.626\pm 0.019$ | $0.635\pm 0.014$ | $0.638\pm 0.015$
site-rshSCP | $0.667\pm 0.016$ | $0.672\pm 0.020$ | $0.673\pm 0.014$ | $0.687\pm 0.015$
comp-rshSCP | $\mathbf{0.703\pm 0.018}$ | $\mathbf{0.727\pm 0.015}$ | $\mathbf{0.728\pm 0.016}$ | $\mathbf{0.731\pm 0.013}$
Table 2: Five fold cross validation for $k_{2}=4$ (mean $\pm$ standard deviation).  Method | $k_{1}=5$ | $k_{1}=10$ | $k_{1}=15$ | $k_{1}=20$
---|---|---|---|---
hSCP | $0.579\pm 0.036$ | $0.583\pm 0.029$ | $0.585\pm 0.014$ | $0.593\pm 0.024$
dem-rshSCP | $0.618\pm 0.038$ | $0.626\pm 0.027$ | $0.635\pm 0.036$ | $0.638\pm 0.032$
site-rshSCP | $0.620\pm 0.037$ | $0.632\pm 0.050$ | $0.634\pm 0.038$ | $0.648\pm 0.026$
comp-rshSCP | $\mathbf{0.651\pm 0.031}$ | $\mathbf{0.680\pm 0.024}$ | $\mathbf{0.681\pm 0.042}$ | $\mathbf{0.689\pm 0.047}$
Table 3: Leave one site accuracy for $k_{2}=4$ (mean $\pm$ standard
deviation).

Table 2 and Table 3 show results of five-fold and leave one site cross-
validation. From the results, we can see the baseline method’s dull prediction
performance, which shows the difficulty of the challenge we are tackling. It
can be seen that dis-rshSCP has a better performance than vanilla hSCP. Out of
demographics and site as a covariate, we see that reducing site variability
has a major impact on the prediction performance compared to reducing
demographics information. The best performance is achieved when we remove both
demographics and site information which suggests that removing heterogeneity
from the data can help improve the predictability power of these components,
thus giving more reliable components discriminative of MDD.

###  5.2 Reproducibility

We have shown that our method can extract components with high predictability
power. This section shows that these components are highly reproducible and
generalizable, which is important for the future application of our framework.
We use split-sample reproducibly to measure the generalizability of the
components, which measures how likely a set of components is replicable across
the same population. First the optimal parameters are selected based on the
highest five fold cross validation, then the split sample reproducibility is
computed by dividing the dataset into two random splits with the same
stratification and calculating the correlation between components derived from
the two splits. Table 4 shows split sample reproducibility of the components
extracted from the different methods. In all the experiment, results are
generated by computing reproducibility over $20$ runs. It can be seen that the
results are similar to prediction performance results, i.e., the components
are highly reproducible, and reducing demographics and site information helps
in improving the reproducibility. This suggests even after adding prediction
loss, the method can find highly reproducible components.

Method | $k_{1}=5$ | $k_{1}=10$ | $k_{1}=15$ | $k_{1}=20$
---|---|---|---|---
hSCP | $0.752\pm 0.032$ | $0.713\pm 0.037$ | $0.688\pm 0.039$ | $0.621\pm 0.021$
dem-rshSCP | $0.781\pm 0.024$ | $0.749\pm 0.045$ | $0.739\pm 0.038$ | $0.688\pm 0.032$
site-rshSCP | $0.782\pm 0.031$ | $0.739\pm 0.032$ | $0.728\pm 0.025$ | $0.661\pm 0.028$
comp-rshSCP | $\mathbf{0.805\pm 0.029}$ | $\mathbf{0.761\pm 0.026}$ | $\mathbf{0.725\pm 0.035}$ | $\mathbf{0.692\pm 0.031}$
Table 4: Split sample reproducibility for $k_{2}=4$ (mean $\pm$ standard
deviation).

(a)

(b)

(c)

(d)

(e)

Figure 2: (1) SMN, DMN and SN, (2) DMN and CEN, (3) VN, (4) DMN and (5) VN and
CEN. We use red and blue colors to represent two different sub-networks which
are anti-correlated with each other and regions with a colored part is
correlated among themselves.

###  5.3 Analysis of Components

We select $k_{1}$ based on optimal tradeoff between classification accuracy
and reproducibility, the result of which are provided in Table 2, 3 and 4. We
can see from the Table 2 and 3 that classification accuracy starts to plateau
at $k_{1}=10$, but if we look at reproducibility, it continues to decrease
linearly. Thus we choose $k_{1}$ to have good classification accuracy without
losing reproducibility of the components. Our aim here is to analyze
components that are most predictive of MDD. For this, the components are
selected on basis of $\mathbf{b}$ in the logistic regression (6). We perform a
hypothesis test with the null hypothesis being that the $k$th component is not
discriminative of MDD, i.e., ${b}_{k}$ is $0$ and $\text{p-value}<0.05$. If a
component has a significant positive value of ${b}_{k}$, then it is likely to
have a higher weight in the MDD population than the healthy population and
vice versa. Each column of $\mathbf{W}$ contains a component storing
information about a set of co-activated regions which can be positively or
negatively correlated with each other, and we map normalized values of each
column onto corresponding regions. On the basis of hypothesis test, we
obtained $5$ significant components out of $10$.

Figure 2 displays components most predictive of MDD.The color bar indicates
the contribution of each region in the components. Regions colored blue are
anti-correlated with displayed areas in red. The extracted components store
information between various parts of the human brain, i.e., whether the
regions are correlated or anti-correlated with each other. These regions could
be smaller parts of the brain or clustered to form resting-state functional
networks. Component $1$ comprises of regions of Somatomotor Network (SMN),
Default Mode Network (DMN) and Salience Network (SN), where regions of DMN and
SMN are anti-correlated with SN. Component $2$ comprises of regions of DMN and
Central Executive Network (CEN). Component $3$ and $4$ consists of regions of
DMN and Visual Network (VN), and component $5$ consists of regions of VN and
CEN. In addition to the $5$ fine scale significant components, we also
recovered a significant hierarchical component comprising of $5$ and $2$ shown
in Figure 3. The strength of each component in MDD and healthy individuals is
shown in Figure 4. We observe decreased representation components $2$, $3$ and
$5$ comprising the DMN, CEN and VN in MDD subjects compared to healthy.
Increased representation is seen in the components comprising SMN, DMN and SN
in MDD subjects compared to healthy.

styleA/.style=top color=white, bottom color=white, styleB/.style=top
color=white, bottom color=red!20, delay=content/.wrap value=##1
, for tree= edge path=[edge] (@parentname.parent anchor) – +(0,-16pt)-|
(name.child anchor)edge label; , forked edges, s sep=20mm [  [5 ] [2 ] ] ]

Figure 3: Hierarchical component having significant predictive power. Here the
coarse-scale component is comprised of component $5$ and $2$ consisting of
regions of DMN, CEN, and VN.

(a) Heterogeneity captured by MDD subjects

(b) Heterogeneity captured by healthy subjects

Figure 4: Heterogeneity captured by components $1-5$. The color represents the
strength ($\mathbf{\Lambda}_{1}$) of each component in each individual. We
normalized $\mathbf{\Lambda}_{1}$ for comparison purposes with $1$ being the
highest value and $0$ lowest.

##  6 Discussion

Our aim was to identify sparse hierarchical connectivity patterns
discriminative of MDD and we achieved this with three coupled loss functions:
1) representation learning loss, 2) discriminative loss, and 3) adversarial
loss. Our model cleverly exploits the rs-fMRI correlation matrix structure to
extract sparse patterns through the representation learning loss. The model
also serves as a dimensionality reduction technique helping to extract low-
rank sparse decomposition. The classification loss in the model forces the
decomposition to extract MDD specific group-level patterns. Notice that there
is a slight tradeoff in classification performance at the expense of
representation learning loss. We highlight this as it is essential for
exploration.

This work provides proof-of-principle analysis focusing on age, sex, and site
as diversity factors in the dataset. Adversarial loss helps reduce these
factors, resulting in improved reproducibility, generalizability and
prediction performance. Our results show that the largest improvement is
acheived when accounting for all factors instead of focusing on just one.
However, other factors could be considered for future work, such as
comorbidity, ethnicity/race, open vs. closed eye during fMRI acquisition, etc.
These factors could help further improve pattern stability and classification
performance.

Here, we used split-sample reproducibility along with cross-validation
accuracy to extract components with the aim to not just have high predictive
power but also high reproducibility. We note here that after a certain number
of components, the addition of more components does not increase the accuracy.
However, the addition of components decreases the reproducibility, which can
be attributed to the model capturing noisy components which are not predictive
of MDD. Future investigations could look at various other reproducibility
measures and classification metrics, and analyze optimality of components in
multiple settings to assess the impact on downstream analysis.

We identified $5$ components that are highly predictive of MDD status. These
components store information about the inter- and intra-connectivity within
networks and suggest that MDD is characterized by disruptions in the following
networks: intra-connectivity in the visual and default mode networks, inter-
connectivity between the visual and central executive networks and inter-
connectivity between the salience, default mode and central executive network.
These discriminative patterns revealed by our framework are consistent with
the recent literature on changes in functional connectivity patterns in MDD.

Impaired visual perception has been found in patients with MDD. It is
considered an important aspect of the disease, whereby there is a positive
correlation between the degree of visual disruption and the severity of
symptoms Song et al. (2021) and reduced visual network connectivity Zeng et
al. (2012), Veer et al. (2010). Studies suggest that perceptual impairments
are linked to abnormal cortical processing and disrupted neurotransmitter
systems whilst retinal processing remains intact Nikolaus et al. (2012),
Salmela et al. (2021).

The DMN has frequently been implicated in MDD pathophysiology due to its role
in producing negative, self-referential, ruminative thoughts Hamilton et al.
(2015). There are previous reports of both hypo- and hyper-connectivity within
the DMN in depression Kaiser et al. (2015), Tozzi et al. (2021) but Liang et
al. (2020) suggest that depression is characterized by two subgroups of
patients exhibiting opposing dysfunctional DMN connectivity. These
inconsistent findings could be due to sample variations in symptom profiles as
variability in connectivity within the DMN is positively correlated with
levels of ruminative thoughts Wise et al. (2017) while hypo-connectivity has
been associated with symptom severity in recurrent MDD Yan et al. (2019).

The inconsistent abnormalities of DMN connectivity in depression suggest that
it could instead be the interplay between the DMN and other networks that
leads to the variety of symptoms observed in depression. Indeed, our patient
sample highly expressed abnormal connectivity across three networks, the CEN,
DMN and SN, which have been put forward as being part of a triple network
model of psychopathology Menon (2011). In this model, aberrant saliency
attribution within the SN weakens the engagement of the CEN and disengagement
of the DMN, leading to cognitive and emotional problems. Therefore,
identifying hierarchical connectivity patterns is crucial to understanding the
interaction between networks that give rise to disease. These findings show
that the proposed method can extract meaningful components with high
reproducibility and clinical relevance without traditional seed-based methods,
which rely too heavily on a priori regions of interest. In a nutshell, these
findings could further our understanding of MDD from a functional network
perspective.

There are several future directions from methodological and clinical
perspectives. First, the model robustness can be improved by introducing
masking of correlation matrices. It has been shown that masking of features
Devlin et al. (2019) while learning can improve the robustness of the model
and its predictability power. Second, instead of reducing age and sex related
heterogeneity, one could disentangle components and learn age, sex and MDD
specific components. Our method is limited to finding effects of MDD on brain
connectivity; in the future, we will use the proposed approach to study the
effects of antidepressant medications on brain connectivity in MDD in resting-
state or task-based fMRI Gudayol-Ferré et al. (2015), Brakowski et al. (2017).
Another important direction would be to combine structural connectivity
information using Diffusion Tensor Imaging (DTI) in our optimization model.
Unifying structural, functional and disease information would give a more
comprehensive view of neurobiological abnormalities and altered brain
functioning and improve MDD diagnosing ability.

Our proposed model has a few weaknesses, adding directions for future work.
First, we only consider logistic regression as our classification model, the
results of which might be sub-optimal. Our method can be modified to include
various classifiers, and a comparison study can be performed to find the
optimal classifier with more emphasis on the classification performances.
Here, using SVM like optimization models might be more straightforward than
incorporating XGBoost like models. In this study, our focus was on the
interpretation of the components; we did not evaluate the reproducibility and
classification accuracy for broad values of $k_{1}$ and $k_{2}$. Future models
could benefit from a thorough investigation of this shortcoming. However, this
analysis is beyond the scope of this work. Second, the robust to covariate
model is a multi-task learning model whose loss weights are manually selected.
More advanced techniques such as the Pareto multi task learning model Lin et
al. (2019) and balanced multi task learning framework Liang and Zhang (2020)
can be used to improve the results. Our study lumps all the MDD patients
together and then analyzes the change in functional networks, but previous
studies have signaled that it is a highly heterogeneous psychiatric disorder
Hyman (2008), Miller (2010). Clustering approaches can be used on
$\mathbf{\Lambda}$ to find subtypes of MDD, the analysis of which is beyond
the scope of the current paper.

##  7 Conclusion

This work presents an effective matrix decomposition strategy to combine rs-
fMRI data with clinical information. Our framework is completely data-
dependent and makes minimal assumptions about the data. We extended the method
of reducing site effects in hSCP by adding additional loss terms for reducing
age and sex effects to estimate robust components discriminative of MDD. We
added a discriminator to extract components that represent the MDD population.
The problem is formulated as a minimax non-convex optimization problem and is
solved using adaptive gradient descent. Experimentally, using a pooled dataset
from five different sites, we showed that reducing heterogeneity introduced by
age, sex, and site could improve the prediction capability of the components,
which is validated using fivefold and leave one site out cross-validation. Our
framework robustly identifies brain patterns characterizing MDD and provides
an understanding of the manifestation of the disorder from a functional
networks perspective. Moreover, our evaluation on a large multi-site dataset
validates the reproducibility and generalizability of the framework. In
addition, our model is not limited to MDD and can be easily adapted to other
disorders such as ASD, ADHD, etc. Moreover, it can easily incorporate other
models outside the medical domain, provided we have access to valid network
measures as an input. This greatly broadens the method’s applicability to
numerous applications from varied fields.

## References

  * Abraham et al. (2017) Abraham, Alexandre; Milham, Michael P; Di Martino, Adriana; Craddock, R Cameron; Samaras, Dimitris; Thirion, Bertrand, and Varoquaux, Gael.  Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example.  _NeuroImage_ , 147:736–745, 2017.
  * Akiki and Abdallah (2019) Akiki, Teddy J and Abdallah, Chadi G.  Determining the hierarchical architecture of the human brain using subject-level clustering of functional networks.  _Scientific reports_ , 9(1):1–15, 2019.
  * Al-Sharoa et al. (2018) Al-Sharoa, Esraa; Al-Khassaweneh, Mahmood, and Aviyente, Selin.  Tensor based temporal and multilayer community detection for studying brain dynamics during resting state fmri.  _IEEE Transactions on Biomedical Engineering_ , 66(3):695–709, 2018.
  * Alexander et al. (2017) Alexander, Lindsay M; Escalera, Jasmine; Ai, Lei; Andreotti, Charissa; Febre, Karina; Mangone, Alexander; Vega-Potler, Natan; Langer, Nicolas; Alexander, Alexis; Kovacs, Meagan, and others, .  An open resource for transdiagnostic research in pediatric mental health and learning disorders.  _Scientific data_ , 4(1):1–26, 2017.
  * Alfaro-Almagro et al. (2021) Alfaro-Almagro, Fidel; McCarthy, Paul; Afyouni, Soroosh; Andersson, Jesper LR; Bastiani, Matteo; Miller, Karla L; Nichols, Thomas E, and Smith, Stephen M.  Confound modelling in uk biobank brain imaging.  _NeuroImage_ , 224:117002, 2021.
  * Arbabshirani et al. (2017) Arbabshirani, Mohammad R; Plis, Sergey; Sui, Jing, and Calhoun, Vince D.  Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls.  _Neuroimage_ , 145:137–165, 2017.
  * Ashourvan et al. (2019) Ashourvan, Arian; Telesford, Qawi K; Verstynen, Timothy; Vettel, Jean M, and Bassett, Danielle S.  Multi-scale detection of hierarchical community architecture in structural and functional brain networks.  _PLoS One_ , 14(5):e0215520, 2019.
  * Batmanghelich et al. (2011) Batmanghelich, Nematollah K; Taskar, Ben, and Davatzikos, Christos.  Generative-discriminative basis learning for medical imaging.  _IEEE transactions on medical imaging_ , 31(1):51–69, 2011.
  * Belmaker and Agam (2008) Belmaker, Robert H and Agam, Galila.  Major depressive disorder.  _New England Journal of Medicine_ , 358(1):55–68, 2008.
  * Benkarim et al. (2021) Benkarim, Oualid; Paquola, Casey; Park, Bo-yong; Kebets, Valeria; Hong, Seok-Jun; de Wael, Reinder Vos; Zhang, Shaoshi; Yeo, BT Thomas; Eickenberg, Michael; Ge, Tian, and others, .  The cost of untracked diversity in brain-imaging prediction.  _bioRxiv_ , 2021.
  * Bhaumik et al. (2017) Bhaumik, Runa; Jenkins, Lisanne M; Gowins, Jennifer R; Jacobs, Rachel H; Barba, Alyssa; Bhaumik, Dulal K, and Langenecker, Scott A.  Multivariate pattern analysis strategies in detection of remitted major depressive disorder using resting state functional connectivity.  _NeuroImage: Clinical_ , 16:390–398, 2017.
  * Biswal et al. (2010) Biswal, Bharat B; Mennes, Maarten; Zuo, Xi-Nian; Gohel, Suril; Kelly, Clare; Smith, Steve M; Beckmann, Christian F; Adelstein, Jonathan S; Buckner, Randy L; Colcombe, Stan, and others, .  Toward discovery science of human brain function.  _Proceedings of the National Academy of Sciences_ , 107(10):4734–4739, 2010.
  * Brakowski et al. (2017) Brakowski, Janis; Spinelli, Simona; Dörig, Nadja; Bosch, Oliver Gero; Manoliu, Andrei; Holtforth, Martin Grosse, and Seifritz, Erich.  Resting state brain network function in major depression–depression symptomatology, antidepressant treatment effects, future research.  _Journal of Psychiatric Research_ , 92:147–159, 2017.
  * Caruana (1997) Caruana, Rich.  Multitask learning.  _Machine learning_ , 28(1):41–75, 1997.
  * Casey et al. (2018) Casey, BJ; Cannonier, Tariq; Conley, May I; Cohen, Alexandra O; Barch, Deanna M; Heitzeg, Mary M; Soules, Mary E; Teslovich, Theresa; Dellarco, Danielle V; Garavan, Hugh, and others, .  The adolescent brain cognitive development (abcd) study: imaging acquisition across 21 sites.  _Developmental cognitive neuroscience_ , 32:43–54, 2018\.
  * Chen et al. (2016) Chen, Heng; Duan, Xujun; Liu, Feng; Lu, Fengmei; Ma, Xujing; Zhang, Youxue; Uddin, Lucina Q, and Chen, Huafu.  Multivariate classification of autism spectrum disorder using frequency-specific resting-state functional connectivity—a multi-center study.  _Progress in Neuro-Psychopharmacology and Biological Psychiatry_ , 64:1–9, 2016.
  * Craddock et al. (2009) Craddock, R Cameron; Holtzheimer III, Paul E; Hu, Xiaoping P, and Mayberg, Helen S.  Disease state prediction from resting state functional connectivity.  _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , 62(6):1619–1628, 2009.
  * Devlin et al. (2019) Devlin, Jacob; Chang, Ming-Wei; Lee, Kenton, and Toutanova, Kristina.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, 2019.
  * Di Martino et al. (2017) Di Martino, Adriana; O’connor, David; Chen, Bosi; Alaerts, Kaat; Anderson, Jeffrey S; Assaf, Michal; Balsters, Joshua H; Baxter, Leslie; Beggiato, Anita; Bernaerts, Sylvie, and others, .  Enhancing studies of the connectome in autism using the autism brain imaging data exchange ii.  _Scientific data_ , 4(1):1–15, 2017.
  * Dong et al. (2019) Dong, Qinglin; Ge, Fangfei; Ning, Qiang; Zhao, Yu; Lv, Jinglei; Huang, Heng; Yuan, Jing; Jiang, Xi; Shen, Dinggang, and Liu, Tianming.  Modeling hierarchical brain networks via volumetric sparse deep belief network.  _IEEE Transactions on Biomedical Engineering_ , 67(6):1739–1748, 2019.
  * Du and Fan (2013) Du, Yuhui and Fan, Yong.  Group information guided ica for fmri data analysis.  _Neuroimage_ , 69:157–197, 2013.
  * Duncan and Northoff (2013) Duncan, Niall W and Northoff, Georg.  Overview of potential procedural and participant-related confounds for neuroimaging of the resting state.  _Journal of psychiatry & neuroscience: JPN_, 38(2):84, 2013.
  * D’Souza et al. (2020) D’Souza, Niharika Shimona; Nebel, Mary Beth; Wymbs, Nicholas; Mostofsky, Stewart H, and Venkataraman, Archana.  A joint network optimization framework to predict clinical severity from resting state functional mri data.  _NeuroImage_ , 206:116314, 2020.
  * Eavani et al. (2013) Eavani, Harini; Satterthwaite, Theodore D; Gur, Raquel E; Gur, Ruben C, and Davatzikos, Christos.  Unsupervised learning of functional network dynamics in resting state fmri.  In _International conference on information processing in medical imaging_ , pages 426–437. Springer, 2013.
  * Eavani et al. (2014) Eavani, Harini; Satterthwaite, Theodore D; Gur, Raquel E; Gur, Ruben C, and Davatzikos, Christos.  Discriminative sparse connectivity patterns for classification of fmri data.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 193–200. Springer, 2014.
  * Eavani et al. (2015) Eavani, Harini; Satterthwaite, Theodore D; Filipovych, Roman; Gur, Raquel E; Gur, Ruben C, and Davatzikos, Christos.  Identifying sparse connectivity patterns in the brain using resting-state fmri.  _Neuroimage_ , 105:286–299, 2015.
  * Ferrarini et al. (2009) Ferrarini, Luca; Veer, Ilya M; Baerends, Evelinda; van Tol, Marie-José; Renken, Remco J; van der Wee, Nic JA; Veltman, Dirk J; Aleman, André; Zitman, Frans G; Penninx, Brenda WJH, and others, .  Hierarchical functional modularity in the resting-state human brain.  _Human brain mapping_ , 30(7):2220–2231, 2009\.
  * First et al. (2004) First, Michael B; France, Allen, and Pincus, Harold Alan.  _DSM-IV-TR guidebook._ American Psychiatric Publishing, Inc., 2004.
  * Frances et al. (1995) Frances, Allen; First, Michael B, and Pincus, Harold Alan.  _DSM-IV guidebook._ American Psychiatric Association, 1995.
  * Godlewska et al. (2014) Godlewska, Beata R; Hasselmann, Helge WW; Igoumenou, Artemis; Norbury, Ray, and Cowen, Philip J.  Short-term escitalopram treatment and hippocampal volume.  _Psychopharmacology_ , 231(23):4579–4581, 2014\.
  * Godlewska et al. (2018) Godlewska, Beata R; Browning, Michael; Norbury, Ray; Igoumenou, Artemis; Cowen, Philip J, and Harmer, Catherine J.  Predicting treatment response in depression: the role of anterior cingulate cortex.  _International Journal of Neuropsychopharmacology_ , 21(11):988–996, 2018.
  * Greicius et al. (2007) Greicius, Michael D; Flores, Benjamin H; Menon, Vinod; Glover, Gary H; Solvason, Hugh B; Kenna, Heather; Reiss, Allan L, and Schatzberg, Alan F.  Resting-state functional connectivity in major depression: abnormally increased contributions from subgenual cingulate cortex and thalamus.  _Biological psychiatry_ , 62(5):429–437, 2007\.
  * Griffanti et al. (2014) Griffanti, Ludovica; Salimi-Khorshidi, Gholamreza; Beckmann, Christian F; Auerbach, Edward J; Douaud, Gwenaëlle; Sexton, Claire E; Zsoldos, Enikő; Ebmeier, Klaus P; Filippini, Nicola; Mackay, Clare E, and others, .  Ica-based artefact removal and accelerated fmri acquisition for improved resting state network imaging.  _Neuroimage_ , 95:232–247, 2014.
  * Gudayol-Ferré et al. (2015) Gudayol-Ferré, Esteve; Peró-Cebollero, Maribel; González-Garrido, Andrés A, and Guàrdia-Olmos, Joan.  Changes in brain connectivity related to the treatment of depression measured through fmri: a systematic review.  _Frontiers in human neuroscience_ , 9:582, 2015.
  * Hamilton et al. (2015) Hamilton, J Paul; Farmer, Madison; Fogelman, Phoebe, and Gotlib, Ian H.  Depressive rumination, the default-mode network, and the dark matter of clinical neuroscience.  _Biological psychiatry_ , 78(4):224–230, 2015\.
  * Hu et al. (2018) Hu, Xintao; Huang, Heng; Peng, Bo; Han, Junwei; Liu, Nian; Lv, Jinglei; Guo, Lei; Guo, Christine, and Liu, Tianming.  Latent source mining in fmri via restricted boltzmann machine.  _Human brain mapping_ , 39(6):2368–2380, 2018\.
  * Huang et al. (2017) Huang, Heng; Hu, Xintao; Zhao, Yu; Makkie, Milad; Dong, Qinglin; Zhao, Shijie; Guo, Lei, and Liu, Tianming.  Modeling task fmri data via deep convolutional autoencoder.  _IEEE transactions on medical imaging_ , 37(7):1551–1561, 2017.
  * Hyman (2008) Hyman, Steven E.  A glimmer of light for neuropsychiatric disorders.  _Nature_ , 455(7215):890, 2008.
  * Iwabuchi et al. (2015) Iwabuchi, Sarina J; Krishnadas, Rajeev; Li, Chunbo; Auer, Dorothee P; Radua, Joaquim, and Palaniyappan, Lena.  Localized connectivity in depression: a meta-analysis of resting state functional imaging studies.  _Neuroscience & Biobehavioral Reviews_, 51:77–86, 2015\.
  * Jenkinson et al. (2012) Jenkinson, Mark; Beckmann, Christian F; Behrens, Timothy EJ; Woolrich, Mark W, and Smith, Stephen M.  Fsl.  _Neuroimage_ , 62(2):782–790, 2012.
  * Jovicich et al. (2016) Jovicich, Jorge; Minati, Ludovico; Marizzoni, Moira; Marchitelli, Rocco; Sala-Llonch, Roser; Bartrés-Faz, David; Arnold, Jennifer; Benninghoff, Jens; Fiedler, Ute; Roccatagliata, Luca, and others, .  Longitudinal reproducibility of default-mode network connectivity in healthy elderly participants: a multicentric resting-state fmri study.  _Neuroimage_ , 124:442–454, 2016.
  * Kaiser et al. (2015) Kaiser, Roselinde H; Andrews-Hanna, Jessica R; Wager, Tor D, and Pizzagalli, Diego A.  Large-scale network dysfunction in major depressive disorder: a meta-analysis of resting-state functional connectivity.  _JAMA psychiatry_ , 72(6):603–611, 2015.
  * Kessler (2012) Kessler, Ronald C.  The costs of depression.  _Psychiatric Clinics_ , 35(1):1–14, 2012.
  * Kostro et al. (2014) Kostro, Daniel; Abdulkadir, Ahmed; Durr, Alexandra; Roos, Raymund; Leavitt, Blair R; Johnson, Hans; Cash, David; Tabrizi, Sarah J; Scahill, Rachael I; Ronneberger, Olaf, and others, .  Correction of inter-scanner and within-subject variance in structural mri based automated diagnosing.  _NeuroImage_ , 98:405–415, 2014.
  * Lee et al. (2010) Lee, Kangjoo; Tak, Sungho, and Ye, Jong Chul.  A data-driven sparse glm for fmri analysis using sparse dictionary learning with mdl criterion.  _IEEE Transactions on Medical Imaging_ , 30(5):1076–1089, 2010.
  * Liang and Zhang (2020) Liang, Sicong and Zhang, Yu.  A simple general approach to balance task difficulty in multi-task learning.  _arXiv preprint arXiv:2002.04792_ , 2020.
  * Liang et al. (2020) Liang, Sugai; Deng, Wei; Li, Xiaojing; Greenshaw, Andrew J; Wang, Qiang; Li, Mingli; Ma, Xiaohong; Bai, Tong-Jian; Bo, Qi-Jing; Cao, Jun, and others, .  Biotypes of major depressive disorder: Neuroimaging evidence from resting-state default mode network patterns.  _NeuroImage: Clinical_ , 28:102514, 2020.
  * Lin et al. (2019) Lin, Xi; Zhen, Hui-Ling; Li, Zhenhua; Zhang, Qing-Fu, and Kwong, Sam.  Pareto multi-task learning.  _Advances in neural information processing systems_ , 32:12060–12070, 2019.
  * Liu et al. (2013) Liu, Feng; Guo, Wenbin; Liu, Ling; Long, Zhiliang; Ma, Chaoqiong; Xue, Zhimin; Wang, Yifeng; Li, Jun; Hu, Maorong; Zhang, Jianwei, and others, .  Abnormal amplitude low-frequency oscillations in medication-naive, first-episode patients with major depressive disorder: a resting-state fmri study.  _Journal of affective disorders_ , 146(3):401–406, 2013.
  * Long et al. (2015) Long, Mingsheng; Cao, Zhangjie; Wang, Jianmin, and Yu, Philip S.  Learning multiple tasks with multilinear relationship networks.  _arXiv preprint arXiv:1506.02117_ , 2015.
  * Menon (2011) Menon, Vinod.  Large-scale brain networks and psychopathology: a unifying triple network model.  _Trends in cognitive sciences_ , 15(10):483–506, 2011.
  * Miller (2010) Miller, Greg.  Beyond dsm: seeking a brain-based classification of mental illness, 2010\.
  * Mueller et al. (1999) Mueller, Timothy I; Leon, Andrew C; Keller, Martin B; Solomon, David A; Endicott, Jean; Coryell, William; Warshaw, Meredith, and Maser, Jack D.  Recurrence after recovery from major depressive disorder during 15 years of observational follow-up.  _American Journal of Psychiatry_ , 156(7):1000–1006, 1999.
  * Mulders et al. (2015) Mulders, Peter C; van Eijndhoven, Philip F; Schene, Aart H; Beckmann, Christian F, and Tendolkar, Indira.  Resting-state functional connectivity in major depressive disorder: a review.  _Neuroscience & Biobehavioral Reviews_, 56:330–344, 2015\.
  * Navrady et al. (2018) Navrady, LB; Wolters, MK; MacIntyre, DJ; Clarke, Toni-Kim; Campbell, AI; Murray, AD; Evans, KL; Seckl, Jonathan; Haley, Christopher; Milburn, Keith, and others, .  Cohort profile: stratifying resilience and depression longitudinally (stradl): a questionnaire follow-up of generation scotland: Scottish family health study (gs: Sfhs).  _International journal of epidemiology_ , 47(1):13–14g, 2018.
  * Nielsen et al. (2013) Nielsen, Jared A; Zielinski, Brandon A; Fletcher, P Thomas; Alexander, Andrew L; Lange, Nicholas; Bigler, Erin D; Lainhart, Janet E, and Anderson, Jeffrey S.  Multisite functional connectivity mri classification of autism: Abide results.  _Frontiers in human neuroscience_ , 7:599, 2013.
  * Nikolaus et al. (2012) Nikolaus, Susanne; Hautzel, Hubertus; Heinzel, Alexander, and Müller, Hans-Wilhelm.  Key players in major and bipolar depression—a retrospective analysis of in vivo imaging studies.  _Behavioural brain research_ , 232(2):358–390, 2012.
  * Noble et al. (2017) Noble, Stephanie; Scheinost, Dustin; Finn, Emily S; Shen, Xilin; Papademetris, Xenophon; McEwen, Sarah C; Bearden, Carrie E; Addington, Jean; Goodyear, Bradley; Cadenhead, Kristin S, and others, .  Multisite reliability of mr-based functional connectivity.  _Neuroimage_ , 146:959–970, 2017.
  * Otte et al. (2016) Otte, Christian; Gold, Stefan M; Penninx, Brenda W; Pariante, Carmine M; Etkin, Amit; Fava, Maurizio; Mohr, David C, and Schatzberg, Alan F.  Major depressive disorder.  _Nature reviews Disease primers_ , 2(1):1–20, 2016.
  * Potluru and Calhoun (2008) Potluru, Vamsi K and Calhoun, Vince D.  Group learning using contrast nmf: Application to functional and structural mri of schizophrenia.  In _2008 IEEE International Symposium on Circuits and Systems_ , pages 1336–1339. IEEE, 2008.
  * Reddi et al. (2019) Reddi, Sashank J; Kale, Satyen, and Kumar, Sanjiv.  On the convergence of adam and beyond.  _arXiv preprint arXiv:1904.09237_ , 2019.
  * Rosa et al. (2015) Rosa, Maria J; Portugal, Liana; Hahn, Tim; Fallgatter, Andreas J; Garrido, Marta I; Shawe-Taylor, John, and Mourao-Miranda, Janaina.  Sparse network-based models for patient classification using fmri.  _Neuroimage_ , 105:493–506, 2015.
  * Ruder (2017) Ruder, Sebastian.  An overview of multi-task learning in deep neural networks.  _arXiv preprint arXiv:1706.05098_ , 2017.
  * Ruder et al. (2019) Ruder, Sebastian; Bingel, Joachim; Augenstein, Isabelle, and Søgaard, Anders.  Latent multi-task architecture learning.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 4822–4829, 2019.
  * Sahoo and Davatzikos (2021a) Sahoo, Dushyant and Davatzikos, Christos.  Extraction of hierarchical functional connectivity components in human brain using adversarial learning.  _arXiv preprint arXiv:2104.10255_ , 2021a.
  * Sahoo and Davatzikos (2021b) Sahoo, Dushyant and Davatzikos, Christos.  Learning robust hierarchical patterns of human brain across many fmri studies.  _arXiv preprint arXiv:2105.06535_ , 2021b.
  * Sahoo et al. (2018) Sahoo, Dushyant; Honnorat, Nicolas, and Davatzikos, Christos.  Gpu accelerated extraction of sparse granger causality patterns.  In _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ , pages 604–607. IEEE, 2018.
  * Sahoo et al. (2019) Sahoo, Dushyant; Honnorat, Nicolas, and Davatzikos, Christos.  Sparse low-dimensional causal modeling for the analysis of brain function.  In _Medical Imaging 2019: Image Processing_ , volume 10949, page 109492R. International Society for Optics and Photonics, 2019.
  * Sahoo et al. (2020) Sahoo, Dushyant; Satterthwaite, Theodore D, and Davatzikos, Christos.  Hierarchical extraction of functional connectivity components in human brain using resting-state fmri.  _IEEE Transactions on Medical Imaging_ , 2020.
  * Salimi-Khorshidi et al. (2014) Salimi-Khorshidi, Gholamreza; Douaud, Gwenaëlle; Beckmann, Christian F; Glasser, Matthew F; Griffanti, Ludovica, and Smith, Stephen M.  Automatic denoising of functional mri data: combining independent component analysis and hierarchical fusion of classifiers.  _Neuroimage_ , 90:449–468, 2014.
  * Salmela et al. (2021) Salmela, Viljami; Socada, Lumikukka; Söderholm, John; Heikkilä, Roope; Lahti, Jari; Ekelund, Jesper, and Isometsä, Erkki.  Reduced visual contrast suppression during major depressive episodes.  _Journal of Psychiatry and Neuroscience_ , 46(2):E222–E231, 2021.
  * Sener and Koltun (2018) Sener, Ozan and Koltun, Vladlen.  Multi-task learning as multi-objective optimization.  _arXiv preprint arXiv:1810.04650_ , 2018.
  * Shinohara et al. (2017) Shinohara, Russell T; Oh, Jiwon; Nair, Govind; Calabresi, Peter A; Davatzikos, Christos; Doshi, Jimit; Henry, Roland G; Kim, Gloria; Linn, Kristin A; Papinutto, Nico, and others, .  Volumetric analysis from a harmonized multisite brain mri study of a single subject with multiple sclerosis.  _American Journal of Neuroradiology_ , 38(8):1501–1509, 2017.
  * Smith et al. (2009) Smith, Stephen M; Fox, Peter T; Miller, Karla L; Glahn, David C; Fox, P Mickle; Mackay, Clare E; Filippini, Nicola; Watkins, Kate E; Toro, Roberto; Laird, Angela R, and others, .  Correspondence of the brain’s functional architecture during activation and rest.  _Proceedings of the national academy of sciences_ , 106(31):13040–13045, 2009.
  * Smith et al. (2014) Smith, Stephen M; Hyvärinen, Aapo; Varoquaux, Gaël; Miller, Karla L, and Beckmann, Christian F.  Group-pca for very large fmri datasets.  _Neuroimage_ , 101:738–749, 2014.
  * Song et al. (2021) Song, Xue Mei; Hu, Xi-Wen; Li, Zhe; Gao, Yuan; Ju, Xuan; Liu, Dong-Yu; Wang, Qian-Nan; Xue, Chuang; Cai, Yong-Chun; Bai, Ruiliang, and others, .  Reduction of higher-order occipital gaba and impaired visual perception in acute major depressive disorder.  _Molecular Psychiatry_ , pages 1–9, 2021.
  * Stolicyn et al. (2020) Stolicyn, Aleks; Harris, Mathew A; Shen, Xueyi; Barbu, Miruna C; Adams, Mark J; Hawkins, Emma L; de Nooij, Laura; Yeung, Hon Wah; Murray, Alison D; Lawrie, Stephen M, and others, .  Automated classification of depression from structural brain measures across two independent community-based cohorts.  _Human brain mapping_ , 41(14):3922–3937, 2020\.
  * Tozzi et al. (2021) Tozzi, Leonardo; Zhang, Xue; Chesnut, Megan; Holt-Gosselin, Bailey; Ramirez, Carolina A, and Williams, Leanne M.  Reduced functional connectivity of default mode network subsystems in depression: meta-analytic evidence and relationship with trait rumination.  _NeuroImage: Clinical_ , 30:102570, 2021.
  * Trivedi et al. (2016) Trivedi, Madhukar H; McGrath, Patrick J; Fava, Maurizio; Parsey, Ramin V; Kurian, Benji T; Phillips, Mary L; Oquendo, Maria A; Bruder, Gerard; Pizzagalli, Diego; Toups, Marisa, and others, .  Establishing moderators and biosignatures of antidepressant response in clinical care (embarc): Rationale and design.  _Journal of psychiatric research_ , 78:11–23, 2016.
  * Veer et al. (2010) Veer, Ilya M; Beckmann, Christian; Van Tol, Marie-Jose; Ferrarini, Luca; Milles, Julien; Veltman, Dick; Aleman, André; Van Buchem, Mark A; Van Der Wee, Nic JA, and Rombouts, Serge AR.  Whole brain resting-state analysis reveals decreased functional connectivity in major depression.  _Frontiers in systems neuroscience_ , 4:41, 2010.
  * Wise et al. (2017) Wise, T; Marwood, L; Perkins, AM; Herane-Vives, A; Joules, R; Lythgoe, DJ; Luh, WM; Williams, SCR; Young, AH; Cleare, AJ, and others, .  Instability of default mode network connectivity in major depression: a two-sample confirmation study.  _Translational psychiatry_ , 7(4):e1105–e1105, 2017.
  * Wu et al. (2011) Wu, Qi-Zhu; Li, Dong-Ming; Kuang, Wei-Hong; Zhang, Ti-Jiang; Lui, Su; Huang, Xiao-Qi; Chan, Raymond CK; Kemp, Graham J, and Gong, Qi-Yong.  Abnormal regional spontaneous neural activity in treatment-refractory depression revealed by resting-state fmri.  _Human brain mapping_ , 32(8):1290–1299, 2011\.
  * Yamashita et al. (2019) Yamashita, Ayumu; Yahata, Noriaki; Itahashi, Takashi; Lisi, Giuseppe; Yamada, Takashi; Ichikawa, Naho; Takamura, Masahiro; Yoshihara, Yujiro; Kunimatsu, Akira; Okada, Naohiro, and others, .  Harmonization of resting-state functional mri data across multiple imaging sites via the separation of site differences into sampling bias and measurement bias.  _PLoS biology_ , 17(4):e3000042, 2019.
  * Yan et al. (2019) Yan, Chao-Gan; Chen, Xiao; Li, Le; Castellanos, Francisco Xavier; Bai, Tong-Jian; Bo, Qi-Jing; Cao, Jun; Chen, Guan-Mao; Chen, Ning-Xuan; Chen, Wei, and others, .  Reduced default mode network functional connectivity in patients with recurrent major depressive disorder.  _Proceedings of the National Academy of Sciences_ , 116(18):9078–9083, 2019.
  * Ye et al. (2015) Ye, Ming; Yang, Tianliang; Qing, Peng; Lei, Xu; Qiu, Jiang, and Liu, Guangyuan.  Changes of functional brain networks in major depressive disorder: a graph theoretical analysis of resting-state fmri.  _PloS one_ , 10(9):e0133775, 2015.
  * Yu et al. (2018) Yu, Meichen; Linn, Kristin A; Cook, Philip A; Phillips, Mary L; McInnis, Melvin; Fava, Maurizio; Trivedi, Madhukar H; Weissman, Myrna M; Shinohara, Russell T, and Sheline, Yvette I.  Statistical harmonization corrects site effects in functional connectivity measurements from multi-site fmri data.  _Human brain mapping_ , 39(11):4213–4227, 2018\.
  * Zeng et al. (2012) Zeng, Ling-Li; Shen, Hui; Liu, Li; Wang, Lubin; Li, Baojuan; Fang, Peng; Zhou, Zongtan; Li, Yaming, and Hu, Dewen.  Identifying major depression using whole-brain functional connectivity: a multivariate pattern analysis.  _Brain_ , 135(5):1498–1507, 2012.
  * Zhang et al. (2020) Zhang, Wei; Zhao, Shijie; Hu, Xintao; Dong, Qinglin; Huang, Heng; Zhang, Shu; Zhao, Yu; Dai, Haixing; Ge, Fangfei; Guo, Lei, and others, .  Hierarchical organization of functional brain networks revealed by hybrid spatiotemporal deep learning.  _Brain connectivity_ , 10(2):72–82, 2020.
  * Zhang and Yang (2017) Zhang, Yu and Yang, Qiang.  A survey on multi-task learning.  _arXiv preprint arXiv:1707.08114_ , 2017.
  * Zhao et al. (2020a) Zhao, Jianlong; Huang, Jinjie; Zhi, Dongmei; Yan, Weizheng; Ma, Xiaohong; Yang, Xiao; Li, Xianbin; Ke, Qing; Jiang, Tianzi; Calhoun, Vince D, and others, .  Functional network connectivity (fnc)-based generative adversarial network (gan) and its applications in classification of mental disorders.  _Journal of neuroscience methods_ , 341:108756, 2020a.
  * Zhao et al. (2020b) Zhao, Youjin; Niu, Running; Lei, Du; Shah, Chandan; Xiao, Yuan; Zhang, Wenjing; Chen, Ziqi; Lui, Su, and Gong, Qiyong.  Aberrant gray matter networks in non-comorbid medication-naive patients with major depressive disorder and those with social anxiety disorder.  _Frontiers in Human Neuroscience_ , 14:172, 2020b.
  * Zhu et al. (2012) Zhu, Xueling; Wang, Xiang; Xiao, Jin; Liao, Jian; Zhong, Mingtian; Wang, Wei, and Yao, Shuqiao.  Evidence of a dissociation pattern in resting-state default mode network connectivity in first-episode, treatment-naive major depression patients.  _Biological psychiatry_ , 71(7):611–617, 2012\.

##  Appendix A Related Work: Extended

Several different methods have been developed to analyze brain organization
using fMRI data. One of the widely used methods include Independent Component
Analysis Smith et al. (2009), Non-Negative Matrix Factorization Potluru and
Calhoun (2008),Sparse Dictionary Learning Lee et al. (2010), Eavani et al.
(2015) and matrix factorization based Granger Causality Sahoo et al. (2018,
2019). These methods are used to estimate interpretable patterns capturing
functional human brain information. Another commonly used method is multi-
scale community detection methods to understand the hierarchical organization
of the human brain Ferrarini et al. (2009), Al-Sharoa et al. (2018), Akiki and
Abdallah (2019), Ashourvan et al. (2019). The approaches mentioned above have
one or more limitations: 1) the assumption that the components are
independent, 2) the inability to capture a subject-specific representation of
the patterns, and 3) removal of negative edge links in the method because the
negative links are treated as repulsions. Recently, several deep learning
based methods Zhang et al. (2020), Huang et al. (2017), Hu et al. (2018),
Zhang et al. (2020), Dong et al. (2019) have been introduced to estimate
functional networks. These methods have shown promising results in terms of
prediction performances. Still, there are one or more disadvantages: 1) not
capturing negative correlation between nodes in a pattern resulting in the
lacking inference of inter-network connectivity, 2) the inability to capture a
subject-specific representation of the patterns, and 3) “black-box” results
due to non-linearity of the deep learning model causing loss of
interpretability.

Hierarchical Sparse Connectivity Patterns (hSCPs) Sahoo and Davatzikos (2021b)
is a dimensionality reduction method overcoming the disadvantages mentioned
above. The method allows us to investigate hierarchically organized co-
activated brain regions and probe individual representations of brain activity
patterns. We consider hSCP as a method that produces a small set of components
consisting of functionally synchronous brain regions that statistically appear
to have relatively more synchronous activation patterns and hence could be
presumed to be parts of underlying brain networks. These components can
capture interactions between various resting-state functional connectivity
networks, for example, Default Mode Network, Salience Network, etc. However,
these networks need not necessarily be present in each individual or subset of
individuals.

##  Appendix B Algorithm

###  B.1 Alternating Minimization

We employ the same alternating minimization technique as used in previous hSCP
papers Sahoo and Davatzikos (2021b, a), Sahoo et al. (2020) for estimating
model parameters. Here, we optimize the objective function 8 for each variable
using adaptive gradient descent (AMSGrad) Reddi et al. (2019) while holding
estimates of other variables as constants. $\beta_{1}$ and $\beta_{2}$ value
in AMSGrad are kept to be $0.9$ and $0.999$. The gradients of each variable
used in gradient descent is defined in the next section. Alternating
minimization procedure to solve equation 8 is described in Algorithm 1. The
algorithm can be modified for solving equation 9 by commenting out the
gradient descent of $\mathbf{b}$. $\mathcal{U}$ and $\mathbf{V}$ are
initialized using equation $10$ in Sahoo and Davatzikos (2021b) and
$\mathop{\rm svd-initialization}\nolimits$ algorithm Sahoo et al.
(2020)[Algorithm 2] is used to initialize $\mathcal{W}$ and $\mathcal{D}$.
Additional constraints in the optimization problem are met by using below
operators:

  * •

$\mathop{\rm proj}\nolimits_{1}(\mathbf{W},\lambda)$\- Projecting each column
of $\mathbf{W}$ into the intersection of $L_{1}$ and $L_{\infty}$ ball

  * •

$\mathop{\rm proj}\nolimits_{2}$\- Projecting a matrix onto $\mathbb{R}_{+}$
space

  * •

$\mathop{\rm proj}\nolimits_{3}$\- Projecting a vector onto $L_{1}$ ball

Algorithm 1 dis-rshSCP

1:Input: Data $\mathcal{C}$, number of connectivity patterns $k_{1}$, …,
$k_{K}$ and sparsity $\lambda_{1}$, …, $\lambda_{K}$ at different level,
hyperparameters $\mu$, $\gamma_{1}$ and $\gamma_{2}$.

2:Initialize $\mathcal{W}$ and $\mathcal{D}$ using $\mathop{\rm svd-
initialization}\nolimits$

3:Initialize $\mathcal{U}$ and $\mathbf{V}$ using $\mathop{\rm site-
initialization}\nolimits$

4:repeat

5: for $r=1$ to $K$ do

6: if Starting criterion is met then

7: $\zeta\leftarrow\mathop{\rm descent}\nolimits(\zeta)$

8: $\mathbf{b}\leftarrow\mathop{\rm descent}\nolimits(\mathbf{b})$

9: if $r==1$ then

10: $\mathbf{W}_{r}\leftarrow\mathop{\rm
proj}\nolimits_{1}(\mathbf{W}_{r},\lambda_{r})$

11: else

12: $\mathbf{W}_{r}\leftarrow\mathop{\rm proj}\nolimits_{2}(\mathbf{W}_{r})$

13: for $n=1,..,N$ do

14: $\mathbf{\Lambda}_{r}^{n}\leftarrow\mathop{\rm
descent}\nolimits(\mathbf{\Lambda}_{r}^{n})$

15: $\mathbf{\Lambda}_{r}^{n}\leftarrow\mathop{\rm
proj}\nolimits_{2}(\mathbf{\Lambda}_{r}^{n})$

16: for $s=1$ to $S$ do

17: $\mathbf{U}^{s}_{r}\leftarrow\mathop{\rm
descent}\nolimits(\mathbf{U}^{s}_{r})$

18: $\mathbf{V}_{r}\leftarrow\mathop{\rm descent}\nolimits(\mathbf{V}_{r})$

19: $\mathbf{V}_{r}\leftarrow\mathop{\rm
proj}\nolimits_{3}(\mathbf{V}_{r},\mu)$

20:until Stopping criterion is reached

21:Output: $\mathcal{W}$, $\mathcal{L}$ and $\mathbf{b}$

###  B.2 Gradient Calculations

This section defines the gradients used in alternating gradient descent
algorithm, some of which are already defined in Sahoo and Davatzikos (2021b).
Let

$\displaystyle\mathbf{\tilde{W}}_{0}$
$\displaystyle=\mathbf{W}_{0}=\mathbf{I}_{P},\qquad\mathbf{Y}_{r}=\prod_{j=0}^{r}\mathbf{W}_{j},\qquad\mathbf{Q}_{m,n}^{r}=(\prod_{j=1}^{m-r}\mathbf{W}_{j})\mathbf{Z}^{s}_{m-r}\mathbf{\Lambda}^{n}_{m-r}(\prod_{j=1}^{m-r}\mathbf{W}_{j})^{\top},$
$\displaystyle\mathbf{X}^{n}_{r}$
$\displaystyle=\mathbf{\Theta}^{n}-\mathbf{U}^{s}_{r}\mathbf{V}_{r},\qquad\mathbf{H}^{n}_{r}=\mathbf{\Theta}^{n}-(\prod_{j=1}^{r}\mathbf{W}_{j})\mathbf{Z}^{s}_{m-r}\mathbf{\Lambda}_{r}^{n}(\prod_{j=1}^{r}\mathbf{W}_{n})^{\top},$

where $n\in\mathcal{I}_{s}$. Let $J$ be the objective function defined in 8
and then gradient of $J$ with respect to $\mathbf{\Lambda}_{r}^{n}$ is:

$\displaystyle\frac{\partial J}{\partial\mathbf{\Lambda}_{r}^{n}}$
$\displaystyle=\frac{\partial
G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})}{\partial\mathbf{\Lambda}_{r}^{n}}-\gamma_{1}\frac{\partial\mathcal{L}_{1}(\zeta,\mathcal{D},\mathbf{y}_{site},\mathbf{y}_{age},\mathbf{y}_{sex})}{\partial\mathbf{\Lambda}_{r}^{n}}+\gamma_{2}\frac{\partial\mathcal{L}_{2}(\mathbf{b},\mathcal{D},\mathbf{y}_{mdd})}{\partial\mathbf{\Lambda}_{r}^{n}}$
$\displaystyle=(-2\mathbf{Y}_{r}^{\top}\mathbf{X}_{r}^{n}\mathbf{Y}_{r}+2\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r}\mathbf{Z}^{s}_{r}\mathbf{\Lambda}_{r}^{n}\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r})\circ\mathbf{Z}_{r}^{s}+\mathbf{F}.$

where is $\mathbf{F}$ is calculated using automatic differentiation tool
provided by matlab. The gradient of $J$ with respect to $\mathbf{W}_{r}$ is:

$\displaystyle\frac{\partial J}{\partial\mathbf{W}_{r}}$
$\displaystyle=\frac{\partial
G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})}{\partial\mathbf{W}_{r}}$
$\displaystyle=\sum_{n=1}^{N}\sum_{j=r}^{K}-4\mathbf{Y}_{r-1}^{\top}\mathbf{X}_{n}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{Q}_{j,n}^{r}+4\mathbf{Y}_{r-1}^{\top}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{Q}_{j,n}^{r}\mathbf{W}_{r}^{\top}\mathbf{Y}_{r-1}^{\top}\mathbf{Y}_{r-1}\mathbf{W}_{r}\mathbf{Q}_{j,n}^{r}.$

The gradient $J$ with respect to $\mathbf{U}^{s}$ and $\mathbf{V}$ is:

$\displaystyle\frac{\partial J}{\partial\mathbf{U}^{s}_{r}}$
$\displaystyle=\left(\sum_{n=\mathcal{I}_{s}}\left(\mathbf{H}^{n}_{r}-\mathbf{U}^{s}_{r}\mathbf{V}_{r}\right)\mathbf{V}^{\top}_{r}\right)\circ\mathbf{I}_{p}$
$\displaystyle\frac{\partial J}{\partial\mathbf{V}_{r}}$
$\displaystyle=\sum_{s=1}^{S}\sum_{n\in\mathcal{I}_{s}}\mathbf{U}^{s}_{r}\left(\mathbf{H}^{n}_{r}-\mathbf{U}^{s}_{r}\mathbf{V}_{r}\right).$

The gradient of $J$ with respect to $\mathbf{K}_{r}^{s}$ is:

$\displaystyle\frac{\partial J}{\partial\mathbf{K}_{r}^{s}}$
$\displaystyle=\frac{\partial
G(\mathcal{C},\mathcal{W},\mathcal{D},\mathcal{U},\mathbf{V},\mathcal{Z})}{\partial\mathbf{K}_{r}^{s}}$
$\displaystyle=\sum_{n=\mathcal{I}_{s}}(-2\mathbf{Y}_{r}^{\top}\mathbf{X}_{r}^{n}\mathbf{Y}_{r}+2\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r}\mathbf{Z}^{s}_{r}\mathbf{\Lambda}_{r}^{n}\mathbf{Y}_{r}^{\top}\mathbf{Y}_{r})\circ\mathbf{\Lambda}_{r}^{n}.$

###  B.3 Convergence results

Reconstruction error is used to empirically validate the convergence of the
Algorithm 1:

$\displaystyle\frac{\sum_{s=1}^{S}\sum_{n\in\mathcal{I}_{s}}\sum_{r=1}^{K}||\mathbf{\Theta}^{i}-(\prod_{j=1}^{r}\mathbf{W}_{j})\mathbf{Z}^{s}\mathbf{\Lambda}_{r}^{i}(\prod_{j=1}^{r}\mathbf{W}_{j})^{\top}-\mathbf{U}^{s}_{r}\mathbf{V}_{r}||_{F}^{2}}{\sum_{n=1}^{N}K||\mathbf{\Theta}^{i}||_{F}^{2}}.$

Figure 5 shows the reconstruction loss, training error and test error of the
algorithm on one of the folds of 5 fold cross validation. We can see from the
result that the algorithm converges after $400$ iterations. The algorithm runs
without the adversarial and discriminative loss for the first $50$ iterations
therefore no accuracy is determined for that duration, and the reconstruction
loss is constantly decreasing. As the adversarial and discriminative losses
are introduced, the reconstruction loss starts to oscillate and then converges
to a sub optimal loss as compared to if there were no additional losses were
introduced. Whereas the training accuracy keeps on getting better but the test
accuracy converges after $200$ iteration. This shows that the algorithm can
overfit, hence necessary cross-validation is important for selecting the
hyperparamters in the loss function 8.

(a) Training and test accuracy

(b) Reconstruction loss

Figure 5: (1) Training and test accuracy and (2) reconstruction loss of dis-
rshSCP algorithm for $k_{1}=10$ and $k_{2}=4$.

##  Appendix C Multi Task Learning

Input LayerShared LayerTask 1Task 2Task 3Task Specific LayerOutput Layer
Figure 6: Multi task learning framework

Figure 6 shows MTL framework used in our problem. The shared layers contain
the following layers: a fully connected layer with $20$ hidden units, dropout
layer with rate $0.2$, ReLU, a fully-connected layer with $10$ hidden units.
Below are the layer detail for each task-

  * •

Age prediction (Task 1)- A fully connected layer with $10$ hidden units, ReLU,
one output unit

  * •

Site prediction (Task 2)- A fully connected layer with $10$ hidden units,
ReLU, a softmax layer

  * •

Sex prediction (Task 3)- A fully connected layer with $10$ hidden units, ReLU,
a softmax layer