# Evaluating the effects of high-throughput structural neuroimaging predictors
on whole-brain functional connectome outcomes via network-based vector-on-
matrix regression

Tong Lu1, Yuan Zhang2, Vince Lyzinski1, Chuan Bi3, Peter Kochunov4, Elliot
Hong4, Shuo Chen3,5
1Department of Mathematics, University of Maryland, College Park
2 Department of Statistics, The Ohio State University
3Maryland Psychiatric Research Center, School of Medicine, University of
Maryland
4Department of Psychiatry and Behavioral Science, University of Texas Health
Science Center
5Division of Biostatistics and Bioinformatics, School of Medicine, University
of Maryland
<EMAIL_ADDRESS>
###### Abstract

The joint analysis of multimodal neuroimaging data is critical in the field of
brain research because it reveals complex interactive relationships between
neurobiological structures and functions. In this study, we focus on
investigating the effects of structural imaging (SI) features, including white
matter micro-structure integrity (WMMI) and cortical thickness, on the whole
brain functional connectome (FC) network. To achieve this goal, we propose a
network-based vector-on-matrix regression model to characterize the FC-SI
association patterns. We have developed a novel multi-level dense bipartite
and clique subgraph extraction method to identify which subsets of spatially
specific SI features intensively influence organized FC sub-networks. The
proposed method can simultaneously identify highly correlated structural-
connectomic association patterns and suppress false positive findings while
handling millions of potential interactions. We apply our method to a
multimodal neuroimaging dataset of 4,242 participants from the UK Biobank to
evaluate the effects of whole-brain WMMI and cortical thickness on the
resting-state FC. The results reveal that the WMMI on corticospinal tracts and
inferior cerebellar peduncle significantly affect functional connections of
sensorimotor, salience, and executive sub-networks with an average correlation
of 0.81 ($p<0.001$).

Keywords: multi-level graph, brain connectome, structural measures, functional
connectivity, dense clique

##  1 Introduction

Neuroimaging data play a fundamental role in deciphering the operations of the
human brain, the most complex organ. These data come in various modalities,
including magnetic resonance imaging (MRI), diffusion tensor imaging (DTI),
and functional MRI (fMRI). Each modality reveals distinct aspects of the
brain’s structure and functionality. For example, MRI provides high-resolution
images of the brain’s structure, offering valuable physical information such
as size, shape, and cortical thickness. DTI assesses the integrity of white
matter microstructures by calculating fractional anisotropy. The fMRI data
capture dynamic blood flow changes in different brain regions to measure
localized neural activity and functional connections.

In statistical analysis, neuroimaging data are commonly represented in two
forms: _vectors_ (e.g., a list of region-wise cortical thickness measures) and
_association matrices_ (e.g., functional connectivity strengths stored in a
weighted adjacency matrix) (Bullmore and Sporns, , 2009; Wig et al., , 2014;
Fornito et al., , 2016; Wang et al., , 2023). Instead of studying brain
structural imaging (SI) and functional connectivity (FC) data separately,
exploring their intricate interplay could significantly deepen our
understanding of the brain, including its development and aging (Smith et al.,
, 2004; Drevets et al., , 2008; Bowman et al., , 2012; Kemmer et al., , 2018).
For example, brain regions connected by white matter tracts with higher
fractional anisotropy are more likely to demonstrate strong FCs, which, in
turn, can influence cognitive processes such as attention, memory, and
decision-making.

There exists little work on the joint analysis of multi-modal neuroimaging
data despite its clear importance, possibly due to the challenge presented by
ultra-high dimensionality and intertwined data structures. In conventional
brain connectome studies, researchers frequently collect $10^{5}$ FC measures
across hundreds of brain regions and up to $10^{4}$ SI measures, resulting in
billions ($10^{9}$) of FC-SI pairs. This not only creates significant
computational demands but also poses challenges for multiple-testing
correction. Traditional correction methods like the false discovery rate (FDR)
and family-wise error rate (FWER) often yield almost no supra-threshold FC-SI
pairs, as demonstrated by extensive simulation studies. Moreover, FC and SI
display data structures indicative of certain connectomic network space and
spatial dependence, respectively. A joint FC-SI analysis needs to incorporate
these intertwined data structures into comprehensive statistical modeling,
thus producing biologically plausible and interpretable results. Specifically,
our goal in this work is to identify an array of SI variables that
intrinsically influences a group of FCs within a brain connectome sub-network,
rather than those randomly distributed across the whole-brain connectome,
referred to as a systematic pattern of associations. These challenges
underscore the necessity of developing a joint analysis method to address the
complexity of multi-modal neuroimaging data.

Recently, advanced statistical methods have been developed to jointly model
two sets of neuroimaging features by leveraging techniques including
regularization, low rank, and projection models (Wang et al., , 2011; Li et
al., , 2012; Zhu et al., , 2014; Kong et al., , 2019). Many of these methods
have been successfully applied to multi-modal imaging data analysis and
yielded interesting findings (Hayden et al., , 2006; Ball et al., , 2017;
Wehrle et al., , 2020; Zhang et al., , 2022). These statistical methods can be
broadly classified into two categories. The first category uses
regularization-based methods (Zhou and Li, , 2014; Zhu et al., , 2017; Wang et
al., , 2020), where a major limitation of these methods is that the sparsely
selected associations fail to take into account the systematic network-level
impacts of SIs on FC networks. The second category employs dimensional
reduction strategies, such as principal component analysis (PCA) (Hotelling, ,
1933; Jolliffe and Cadima, , 2016; Chachlakis et al., , 2019), which first
projects both FCs and SIs into a handful of top principal components and then
performs regression analysis on these selected components. However, as an
unsupervised dimension reduction technique, PCA-based analysis often extracts
fewer associated principal components of outcomes and predictors, thereby
missing the truly associated FC-SI pairs. Sparse canonical correlation
analysis (sCCA) methods can be considered as an integration of these two
categories and have been widely used in neuroimaging studies (Witten et al., ,
2009; Lin et al., , 2013; Uurtio et al., , 2019). Yet, sCCA methods usually
focus on vector-to-vector association analysis, which may also overlook the
systematic vector-to-network association patterns that are of particular
interest in this work (i.e., the associations between the SI vector and FC
sub-connectome represented as a matrix). To bridge the methodological gap in
modeling vector-to-matrix associations and incorporating latent network
structures, we propose a new multi-level network association method (MOAT) to
systematically investigate the FC-SI association patterns.

Figure 1:  The detection pipeline of systematic FC-SI association patterns by
MOAT. (a) shows three sets of input measures: FC data, SI data, and
confounding variables. (b) shows a heatmap of an association matrix
$\boldsymbol{A}$, with each element quantifying the association between SI and
FC measures (e.g., correlation, test statistics, $-log(p)$ value). A hotter
point indicates a stronger FC-SI association. (c) shows the systematic FC-SI
association patterns identified by MOAT: the first level (c1) depicts the
revealed association patterns between SIs and FCs, where strong associations
are mostly concentrated in two dense sub-networks; (c2) shows the FCs in brain
connectome space from the detected clusters $B_{1}$ and $B_{2}$ in (c1); (c3)
shows the second level, which depicts the organized structure of a graph
formed by the selected FC clusters in (c2). Next, our proposed network-based
inference test is performed on each FC-SI sub-network identified by MOAT. (d)
and (e) shows the 3D visualization of the significant FC-SI network-level
associations that pass the inference test.

Figure 1 presents an overview of the MOAT method, which is constructed based
on a multi-level graph for structural-functional neuroimaging data. The first
level is a bipartite graph that depicts the association patterns between the
SI vector as predictors and vectorized FC outcomes, adjusted for other
confounding covariates (Figure 1(a)). Meanwhile, the second level is a
complete unipartite graph that reconstructs the vectorized FCs back to a
whole-brain connectome network. This multi-level structure enables the
identification of subsets of SIs that systematically impact FC sub-networks.
We have further developed computationally efficient algorithms to extract the
multi-level sub-networks from the full graphs and have proposed a tailored
network-based inference frame to individually test each sub-network with
multiple corrections based on permutation tests. Our method is also compatible
with the existing methods aforementioned (e.g., PCA, CCA). For example,
applying CCA to FCs and SIs in an extracted multi-level sub-network provides
an estimate of association in the context of multiple regressions.

The contributions of this article are three-fold. First, we introduce MOAT, a
novel method that can handle matrix-variate outcomes and vector-variate
predictors. Compared to the existing models for multivariate outcomes and
multivariate predictors (Zhuang et al., , 2017; Wu et al., , 2021; Mihalik et
al., , 2022; Lu et al., , 2023), MOAT can further account for the network
structure within the matrix outcomes and between the outcome-predictor
association patterns. MOAT naturally prohibits most false positive
associations because these associations are more likely distributed sparsely
rather than gathered in organized sub-networks. Secondly, we develop new
algorithms to extract those multi-level sub-networks. The computational load
is low because we developed a tailored greedy peeling algorithm with
multilinear complexity, making our approach compatible with the commonly used
permutation tests that are often computationally intensive. Lastly, we
proposed a novel network-level inference framework, where we utilize novel
test statistics derived based on the multi-level dense subgraph properties in
terms of size and density. This inference framework leads to a simultaneous
enhancement of both sensitivity and specificity by leveraging graph
combinatorial theories.

The rest of this paper is organized as follows. In Section 2, we formally
define the multi-level network structure and present how MOAT works in network
extraction with the network-based inference method. In Section 3, we perform
extensive simulation analyses for method validation and comparison. In Section
4, we apply MOAT to a real structure-function neuroimaging dataset from the UK
Biobank with 4,242 participants to systematically investigate the FC-SI
associations. We conclude with discussions in Section 5.

##  2 Our method

###  2.1 Data structure and problem set up

We collect structural-functional neuroimaging data from independent subjects,
indexed as $\\{1,\ldots,D\\}$. For each subject $d:1\leq d\leq D$, we observe
three sets of measurements:

  1. (i)

Independent variables: a vector of $m$ SI measures
$\boldsymbol{X}^{(d)}=\big{(}x_{1}^{(d)},\ldots,x_{m}^{(d)}\big{)}^{T}$. This
vector characterizes anatomical structures of the brain, such as white matter
microstructure integrity measured by fractional anisotropy from DTI (Mori et
al., , 2008) and region-wise cortical thickness obtained from MRI (Tustison et
al., , 2014).

  2. (ii)

Outcome variables: an adjacency matrix that stores pairwise FC measures
$\boldsymbol{Y}^{(d)}_{n\times n}$ between $n$ brain regions. Each element
$y_{ij}^{(d)}$ of $\boldsymbol{Y}^{(d)}$ represents the strength of functional
connection between brain regions $i$ and $j$ of subject $d$, calculated from
functional imaging data such as resting state fMRI. Thanks to the Brainetcome
Atlas (Fan et al., , 2016), researchers can align the FC brain region
partitions across different participants, thus conveniently, their
$\boldsymbol{Y}$ share a common node set. We model $\boldsymbol{Y}^{(d)}$ as
the outcome variable due to the widely accepted view in neurology that brain
structure determines neural functions (Buckner et al., , 2008; Bai et al., ,
2009; Honey et al., , 2010).

  3. (iii)

Confounding variables:
$\boldsymbol{\eta}^{(d)}=\big{(}\eta_{1}^{(d)},\ldots,\eta_{p}^{(d)}\big{)}^{T}$.
These variables include profiling information such as age, sex, genetics, and
environment that may potentially affect brain functional connectome in
complicated ways.

###  2.2 Multi-level graph representation

We explore the brain structural-functional relationship by considering the
following regression model: for each subject $d\in[1,D]$,

$\displaystyle g(y_{ij}^{(d)})=$ $\displaystyle\leavevmode\nobreak\
\theta_{ij}^{0}+\sum_{k=1}^{m}\beta_{(ij),k}x_{k}^{(d)}+\sum_{p=1}^{P}\alpha_{ij}^{p}\eta_{p}^{(d)},$
(1)

where $g(\cdot)$ is a link function, $\theta_{ij}^{0}$ is the intercept,
$\beta_{(ij),k}$ is the coefficient of the SI measure $x_{k}$, and
$\alpha_{ij}^{p}$ is the coefficient of the nuisance covariate $\eta_{p}$
(Zhang et al., , 2023). The focal parameter of interest in the above
regression model (1) is $\beta_{(ij),k}$, where a nonzero coefficient
$\beta_{(ij),k}\neq 0$ signifies an association between an SI measure $x_{k}$
and the functional connection $y_{ij}$ between brain regions $i$ and $j$.
Consequently, learning the set $\\{\beta_{(ij),k}\neq 0\\}$ allows for the
unveiling of brain-wide association patterns between SIs and FCs.

A multi-level graph model targeting $\\{\beta_{(ij),k}\neq 0\\}$ associations.
To facilitate downstream analysis, we let a matrix
$\boldsymbol{\beta}=\\{\beta_{(ij),k}\\}_{\forall
ijk}\in{\mathbb{R}^{\binom{n}{2}\times m}}$ to denote all SI-FC pair-wise
associations. (Zhang and Xia, , 2018). We build the multi-level graph model
based on the $\binom{n}{2}\times m$ matrix $\boldsymbol{\beta}$. Specifically,
at the first level, we define a bipartite graph $B=(S,F;H)$ to represent the
matrix $\boldsymbol{\beta}$, where $S=\\{1,\ldots,m\\}$ (i.e., $|S|=m$)
constitutes the node set of SI measures; $F=\\{1,\ldots,\binom{n}{2}\\}$
(i.e., $|F|=\binom{n}{2}$) constitutes the node set of FC measures; and $H$
denotes the edge set. Each element $h_{(ij),k}\in H$ signifies a non-zero
association between FC and SI (i.e., $\beta_{(ij),k}\neq 0$). We demonstrate
the first level bipartite graph $B$ in the left panel of Figure 2 . The second
level of the multi-level graph model is a classic graph model reflecting the
whole-brain connectome network, denoted as $G=(V;F)$, where $V$ is the node
set of brain regions with size $|V|=n$, while $F$ is the edge set connecting
brain regions with size $|F|\leq(n-1)/2$. Noticeably, each node $(i,j)$ in $F$
can also be interpreted as an edge in the brain functional connetome network
$G(V;F)$. Thus, $F$ denotes both (i) the node set of the bipartite graph
$B=(S,F;H)$ with $F=\\{f_{(ij)}\\}$, where $f_{(ij)}$ represents a node for
the outcome $Y_{ij}$; (ii) the edge set of $G(V;F)$ with $F=\\{f_{i,j}\\}$,
where $f_{i,j}=1$ indicates that brain areas $i$ and $j$ are connected.

Figure 2:  A toy example of a multi-level graph denoted as $B$ and $G$, where
there exists one sub-network $B_{1}\in B$ exhibiting a significant association
between FC and SI. In Level 1, $B=(S,F;H)$ is an input bipartite graph
consisting of two node sets $S$ and $F$, representing the SI and FC measures
respectively. The edge set $H$ is binary and determined by the regression
coefficient $\boldsymbol{\beta}=\\{\beta_{(ij),k}\\}$ which infers the
association between $y_{ij}$ (FC strength between brain regions $i$ and $j$)
and $x_{k}$ ($k-$th SI measure). In this toy example, there exists a sub-
network $B_{1}=(S_{1},F_{1};H_{1})\subset B$ which contains SI-related edges
$\\{h_{(ij),k}|i,j,k\in B_{1}\\}$ with a much higher density than the rest
$\\{h_{(ij),k}|i,j,k\notin B_{1}\\}$. Given that the node set $F$ represents
the whole-brain functional connectome, each node $f_{(ij)}\in F$ comes from an
edge in a brain connectome graph $G=\\{V;F\\}$ (Level 2 graph), where $V$ is
the node set of brain regions, and $F$ is the edge set of interregional
connections. $\\{f_{(ij)}\\}\in B_{1}$ are not isolated in the brain
connectome; instead, they constitute an organized sub-graph $G_{1}\in G$.

In light of the highly organized brain structures and functions, it is
neurobiologically sensible to model $\\{\beta_{(ij),k}\neq 0\\}$ in organized
association patterns (Craddock et al., , 2013; Bahrami et al., , 2019).
Specifically, we consider that a subset of brain structural predictors jointly
influences connectome outcome variables within a functional subnetwork, which
characterizes a plausible brain structure-function interaction (Zalesky et
al., , 2010; Cao et al., , 2014). Built upon this latent relationship pattern,
we specify that $\\{\beta_{(ij),k}\neq 0\\}$ predominantly concentrates within
specific subgraphs denoted as $\\{B_{c}\\}$ where $c=1,\cdots,C$ and
$B_{c}\subset B$. For simplicity, we are going to illustrate the case where
$c=1$ below.

We specify $B_{1}(S_{1},F_{1})$ as a doubly-dense multi-level subgraph (see
dense graph studied in Tong MOAT, Craddock et al.,  (2013); Wu et al.,
(2021)). At the first level, a subset of SI predictors of $\\{X_{k},k\in
S_{1}\\}$ condensely affect $\\{Y_{ij},ij\in F_{1}\\}$:

$\operatorname{Pr}\left(\beta_{(ij),k}\neq 0\mid(ij)\in F_{1},k\in
S_{1}\right)\gg\operatorname{Pr}\left(\beta_{(ij),k}\neq 0\mid(ij)\notin
F_{1}\textrm{ or }k\notin S_{1}\right).$ (2)

At the second level, a connectomic subnetwork $G_{1}=(V_{1};F_{1})$ is an
edge-induced sub-clique, where the edge subset of interest is
$\\{(i,j):(ij)\in F_{1}\\}$. $G_{1}$ is also dense reflecting that SIs of
$S_{1}$ are associated with connectomic edges in a network rather than
sparsely and randomly distributed in the whole-brain connectome (ref Tong
2023). In Figure 2, we demonstrate a doubly-dense multi-level subgraph $B_{1}$
with red-bold edges. Provided with $\\{B_{c}\\}$, we can express the overall
multi-level graph as (3) and (4) as follows:

$\displaystyle\text{Level 1 (bipartite) network: }B=\cup_{c=1}^{C}B_{c}\cup
B_{0},\leavevmode\nobreak\ B_{c}=(S_{c},F_{c};H_{c});$ (3)
$\displaystyle\text{Level 2 (unipartite) network: }G=\cup_{c=1}^{C}G_{c}\cup
G_{0},\leavevmode\nobreak\ G_{c}=(V_{c};F_{c}),$ (4)

where $B_{c}$ and $G_{c}$ are dense subgraphs in $B$ and $G$ respectively, and
$B_{0}$ and $G_{0}$ are the remaining graphs. Each node set $F_{c}$ in $B_{c}$
corresponds to a subset of edges in the functional connectome $G$, which
induces one or multiple cliques in $G$. For simplicity, we use $G_{c}$ to
denote the clique(s) for the corresponding $B_{c}$. If $B$ is a random graph,
then $B=B_{0}$ and for all $c=1,...,C$, $B_{c}=\emptyset$. Similarly, if $G$
is a random graph, then $G_{c}=\emptyset$. Otherwise, $G_{c}$ represents a
connectome sub-network. In Figure 2, we demonstrate a graphical example of the
multi-level structure of $B_{c}$ and $G_{c}$ when $c=1$. In summary, our
multi-level network model assigns a small proportion of $\\{\beta_{(ij),k}\neq
0\\}$ to structured subnetworks reflecting systematic FC-SI association
patterns. The patterns may not be captured by neither shrinkage regression
models nor clustering/biclustering methods.

###  2.3 Multi-level subnetwork estimation

In practice, neither $\\{\beta_{(ij),k}\neq 0\\}$ nor $\\{B_{c}\\}$ is known
and it is challenging to simultaneously handle billions of FC-SI associations
and estimate $\\{B_{c}\\}$ ($\\{\beta_{(ij),k}\neq 0\\}$) in one big model
such as (1) (Woo et al., , 2014; Mbatchou et al., , 2021; Marek et al., ,
2022). To alleviate the computation burden, we take a divide-and-conquer
approach and run one regression for each $k$, recognizing that both $\theta$
and $\alpha$ may also be different for each $k$. This strategy is commonly
used in large-scale imaging and genetics data analysis (Zalesky et al., ,
2010; Schaid et al., , 2018; Chen et al., , 2023).

Next, we extract the desired dense subgraphs $\\{B_{c}\\}$ based on
$\boldsymbol{X}^{(d)},\boldsymbol{Y}^{(d)},\boldsymbol{\eta}^{(d)}$. Since
$\\{\beta_{(ij),k}\neq 0\\}$ are unknown, we compute an inference measure
$a_{(ij),k}$ as a surrogate to $\beta_{(ij),k}$: each $a_{(ij),k}$ is produced
by the statistical inference of a regression model for $x_{k}$ and $y_{ij}$.
For example, $a_{(ij),k}$ can be the $-\log(p)$ for $\beta_{(ij),k}$, where
$-\log(p)$ is a widely used metric in high-dimensional data analysis, such as
Genome-wide association studies (GWAS) and neuroimaging analysis (Lasky-Su et
al., , 2008; Tang et al., , 2016; Sun et al., , 2022). Now we propose the
following criterion for selecting $(S_{1},F_{1})$:

$\underset{S_{c}\subseteq S,F_{c}\subseteq F,V_{c}\subseteq
V}{\arg\max}\sum_{c=1}^{C}\frac{\sum_{k\in S_{c},(ij)\in
F_{c}}a_{(ij),k}}{\left(\left|S_{c}\right|\binom{|V_{c}|}{2}\right)^{\lambda_{1}/2}}+\frac{\sum_{i,j\in
V_{c},i<j}f_{ij}}{\left|V_{c}\right|^{\lambda_{2}}},$ (5)

where $\lambda_{1},\lambda_{2}$ tune the impacts of the densities of $B_{c}$
and $G_{c}$, respectively. For example, when $\lambda_{1}=2$, the first term
becomes the familiar quantity of subgraph density in network analysis. We
typically search $\lambda_{1}$ within the range of $(1,2)$. Empirically,
setting $\lambda_{1}=2$ usually forces $B_{c}$’s into singletons; while
setting $\lambda_{1}$ below 1 often leads to sparse $B_{c}$’s. Likewise, we
explore the parameter $\lambda_{2}$ within the same interval $(1,2)$.
Deviating from this range for $\lambda_{2}$, either higher or lower, will
yield results similar to those observed for $\lambda_{1}$. Here, we follow the
convention in neuroimaging analysis and select $\lambda_{1}$ and $\lambda_{2}$
using Kullback–Leibler (KL) divergence (Johnson and Sinanovic, , 2001; Yohai,
, 2008; Zhao et al., , 2023). Detailed selection procedures are provided in
Appendix A..

Directly solving (5) requires combinatorial computation. Therefore, we propose
a greedy peeling algorithm as a fast approximation. Our algorithm extends the
greedy algorithm for single-level bipartite subgraphs extraction in Wu et al.,
(2021) and Chekuri et al.,  (2022). We present a condensed version as
Algorithm 1 below, and relegate the detailed step-by-step algorithm to
Appendix B. For each multi-level subgraph $B_{c}$ and $G_{c}$, Algorithm 1
first initializes node sets $S_{c}$ and $F_{c}$ with the nodes $S$ and $F$
from the original full graphs, respectively. It then iteratively removes nodes
with the smallest degree (say, $\tau\in S_{c}$ and $\phi\in F_{c}$) from
either $S$ or $F$ (see Line 5 of the algorithm). At the end of each iteration
$q$, the updated node set $F_{c}^{(q)}$ is used to construct the “level 2”
subgraph $G^{(q)}_{c}(V_{c};F_{c})$, and the corresponding output value of
objective function (5) is recorded. This process of node removal and the
construction of “level 2” graph continues iteratively until all nodes have
been excluded from $S_{c}$ or $F_{c}$, with the termination determined by
whichever node subset is exhausted first. Ultimately, the algorithm returns
the dense subgraph $B_{c}$ that maximizes (5) among all $\\{B_{c}^{(q)}\\}$
(see Line 14 of the algorithm).

The computational complexity of Algorithm 1 is
$\mathcal{O}\Big{(}\mathcal{C}|V|(|S|+|F|)\Big{)}$, where $\mathcal{C}$
depends on the number of the grid search, $|V|=n$, $|S|=m$ and
$|F|=\binom{n}{2}$ are the numbers of regions, SI measures, and FC measures,
respectively. Additionally, Theorem 1 confirms the consistency of multi-level
subgraph detection. In essence, the solution to the objective function (5)
gives a consistent estimation of the true multi-level sub-network structure
represented by ${B_{c}}$ (the set of edge-induced sub-networks). As the sample
size $D\to\infty$, the likelihood of an incorrect edge assignment for
${B_{c}}$ approaches zero.

Algorithm 1 Optimization of objective function (5) \- Condensed pseudo code

1:Input: $B=(S,F;H)$, $G=(V,F)$, $\lambda_{1}$, $\lambda_{2}$

2:output: $\\{B_{c}\\}$

3:Define function [$B_{c}\leavevmode\nobreak\ $ density$(B_{c})$]=
greedy_peeling_MOAT [$B\leavevmode\nobreak\ G\leavevmode\nobreak\
\lambda_{1}\leavevmode\nobreak\ \lambda_{2}]$:

4:Initialize ”level 1” subgraph nodes $S_{c}^{1}\leftarrow
S,F_{c}^{1}\leftarrow F$

5:for $q=1,2,\ldots,m+\binom{n}{2}-1$ do

6: Let $\tau\in S_{c}^{(q)}$ and $\phi\in F_{c}^{(q)}$ be the nodes with
smallest degree

7: If$\sqrt{d}\tau\leq\frac{1}{\sqrt{d}}\phi$, then remove $\tau$ from
$S_{c}^{(q)}$; otherwise, remove $\phi$ from $F_{c}^{(q)}$

8: Next, construct “level 2” subgraph $G_{c}(V_{c};F_{c})$, based on current
$F_{c}$

9: Initialize ”level 2” subgraph nodes $V_{c}^{1}\leftarrow V$

10: for $p=1,2,\ldots,n-1$ do

11: remove node with smallest degree from $V_{c}^{(p)}$, store each
$V_{c}^{(p)}$ and the corresponding $F_{c}^{(p)}$

12: end for

13: Output $G^{(q)}_{c}(V_{c};F_{c})$ that maximizes $\frac{\sum_{i,j\in
V_{c},i<j}f_{ij}}{\left|V_{c}\right|^{\lambda_{2}}}$ among
$\\{V_{c}^{(p)},F_{c}^{(p)}\\}_{p=1}^{n-1}$

14: Replace $F_{c}^{(q)}$ with $F_{c}^{(p)}$

15:end for

16:Output $B_{c}$ that maximizes (5) among
$\\{B_{c}^{(q)}\\}_{q=1}^{m+\binom{n}{2}-1}$

17:while density$(B_{c})>$ density($B$) do

18: Fill median of $\\{a_{(ij),k}|k\in S,(ij)\in F\\}$ into
$\\{a_{(ij),k}|k\in S_{c},(ij)\in F_{c}\\}$, obtain updated $B^{*}$ and
$G^{*}$

19: [$B_{c}\leavevmode\nobreak\ $ density$(B_{c})$]= greedy_peeling_MOAT
[$B^{*}\leavevmode\nobreak\ G^{*}\leavevmode\nobreak\
\lambda_{1}\leavevmode\nobreak\ \lambda_{2}\leavevmode\nobreak\ ]$

20:end while

21:Output all $\\{B_{c}\\}$

######  Theorem 1.

(Consistency of subgraph detection). let
$\mathbf{U}^{*}\in\mathbb{R}^{|S|\times|F|}$ be a matrix storing the true edge
membership in $B$, where each element $u_{(ij),k}^{*}=1$ if
$\beta_{(ij),k}\neq 0$, and $u_{(ij),k}^{*}=0$ otherwise. Similarly, let
$\hat{\mathbf{U}}\in\mathbb{R}^{|S|\times|F|}$ store the edge membership
estimated by optimizing (5), where each element $\hat{u}_{(ij),k}=1$ if
$x_{k}\in\hat{S}_{c}$ and $y_{ij}\in\hat{F}_{c}$; $\hat{u}_{(ij),k}=0$
otherwise. Then, for an arbitrarily small $\epsilon$, when the sample size
$D\to\infty$, we have

$\displaystyle\mathbb{P}(||\mathbf{U}^{*}-\hat{\mathbf{U}}||_{F}<\epsilon)\to
1,$

where $||.||_{F}$ denotes the frobenius norm.

The proof of Theorem 1 is provided in Appendix B.2.

###  2.4 Reduced false positive findings by $B_{c}$

Compared to methods that individually select $\\{\beta_{(ij),k}\\}$ such as
multiple testing approaches, our method selects nonzero $\beta$’s via dense
FC-SI associated sub-network $B_{c}$, which can drastically reduce false
positive findings. Let $\\{\hat{\beta}_{(ij),k}\\}$ denote the set of
estimated association parameters from a sample; then
$\\{\hat{\beta}_{(ij),k}\neq 0|\beta_{(ij),k}=0\\}$ indicate false positive
findings. Following the common practice in neuroimaging and neurobiology
(Margulis and Sagan, , 2000), we assume that false positive associations are
randomly distributed in the brain space. The conventional approach using
individual inference on $\\{\hat{\beta}_{(ij),k}\\}$ may likely select many
false positives $\hat{\beta}_{(ij),k}\neq 0|\beta_{(ij),k}=0$. In contrast,
our method returns few false positives. The reason behind is demonstrated in
the following lemma, which says false positives very rarely form dense
subgraphs of moderate sizes.

######  Lemma 1.

Assume that $B_{c}$ is observed from a random multi-level binary graph with a
bipartite graph $B(S,F;H)$ in Level 1 and a unipartite graph $G(V;F)$ in Level
2. Suppose that $B_{c}$ is a multi-level subgraph that has: (1) Edge density
in $B_{c}$ with $\frac{\sum_{k\in S_{c},(ij)\in
F_{c}}I(\hat{\beta}_{(ij),k}\neq
0|\beta_{(ij),k}=0)}{\left|S_{c}\right|\left|F_{c}\right|}\geq{\gamma_{1}}\in\left(p_{1},1\right)$,
where $p_{1}=\frac{\sum_{k\in S,(ij)\in F}I(\hat{\beta}_{(ij),k}\neq
0|\beta_{(ij),k}=0)}{|S||F|}$ is the proportion of false positive associations
in $B$; (2) Edge density in $G_{c}$ with $\frac{\sum_{i,j\in
V_{c},i<j}I(\hat{\beta}_{(ij),k}\neq
0|\beta_{(ij),k}=0)}{\binom{\left|V_{c}\right|}{2}}\geq{\gamma_{2}}\in\left(p_{2},1\right)$,
where $p_{2}=\frac{\sum_{i,j\in V,i<j}I(\hat{\beta}_{(ij),k}\neq
0|\beta_{(ij),k}=0)}{\binom{|V|}{2}}$ is the proportion of false positive
associations in $G$. Furthermore, let
$m_{0},n_{0}=\Omega\left(\max\left\\{m^{\epsilon},n^{\epsilon}\right\\}\right)$
for some $0<\epsilon<1$, where $\Omega$ denotes a loose lower bound. Then for
sufficiently large $m,n$ with $\zeta(\gamma_{1},p_{1})m_{0}\geq 4\log n(n-1)$,
$\zeta(\gamma_{1},p_{1})n_{0}(n_{0}-1)\geq 16\log m$, and
$\zeta(\gamma_{2},p_{2})n_{0}\geq 4\log n$, we have

$\displaystyle\mathbb{P}\left(\left|S_{c}\right|\geq
m_{0},\left|F_{c}\right|\geq\binom{n_{0}}{2},\left|V_{c}\right|\geq
n_{0}\right)$ $\displaystyle\leq
2mn^{2}(n-1)\cdot\exp\left(-\frac{1}{8}\zeta\left(\gamma_{1},p_{1}\right)m_{0}n_{0}(n_{0}-1)-\frac{1}{4}\zeta\left(\gamma_{2},p_{2}\right)n_{0}^{2}\right),$
(6)

where $\zeta(a,b)=\left\\{\frac{1}{(a-b)^{2}}+\frac{1}{3(a-b)}\right\\}^{-1}$.

Lemma 1 is proved in Appendix B. Lemma 1 states that the probability of
identifying a multi-level subgraph of $\hat{B}_{c}$ composed of false positive
associations $\\{\hat{\beta}_{(ij),k}\neq 0|\beta_{(ij),k}=0\\}$ exponentially
converges to 0 as the sizes and densities of the multi-level sub-network
increase. In practice, the probability of a false positive network with
reasonable size (e.g., $|S|\times|F|=10\times 10$) and sound densities is less
than $10^{-16}$. It is very unlikely that false positive FC-SI associations
$\\{\beta_{(ij),k}\neq 0\\}$ would form a large and dense subgraph $B_{c}$.

###  2.5 Inference for extracted $\hat{B}_{c}$

Recall from Section 2.3 that our study aims to identify specific subsets of
SIs and FCs that exhibit systematic association patterns encoded by $B_{c}$.
Performing Algorithm1 returns a collection of such subgraphs $\hat{B}_{c}$.
Our next goal is to conduct a network-level statistical inference to gauge the
significance of each $\hat{B}_{c}$ with multiple corrections (Goeman et al., ,
2022; Zhang et al., , 2023; Chen et al., , 2023). Roughly speaking, we assess
the statistical significance of each $B_{c}$ by testing:

$\displaystyle\mathbb{H}_{0}$ $\displaystyle:\textrm{$B_{c}$ is not a dense
multi-level subgraph constituted by associated SI-FC pairs;}$
$\displaystyle\mathbb{H}_{a}$ $\displaystyle:\textrm{$B_{c}$ is a dense multi-
level subgraph reflecting systematic FC-SI associations.}$ (7)

More precisely, under $\mathbb{H}_{0}$, the edges of $B_{c}$ is randomly
distributed among all possible pairs. Per Lemma 1, it is rare to observe large
and dense multi-level subgraph $B_{c}$ under the null. Therefore, we can
straightforwardly perform the commonly used permutation testing strategy in
neuroimaging statistics to assess the significance of $\hat{B}_{c}$ while
controlling the FWER (Zalesky et al.,  (2010); Nichols,  (2012); Woo et al.,
(2014). However, our testing object is a multi-level subgraph, which is
different from the voxel-based “clusters” commonly encountered in conventional
cluster-extent inference because the rareness of $\hat{B}_{c}$ is jointly
determined by both the densities and sizes of the dense bipartite and clique
of $\hat{B}_{c}$ instead of a measure of cluster-extent (e.g., the number of
voxels). To address this challenge, we propose a novel test statistic
$\mathcal{T}(\hat{B}_{c})$, proportional to the upper bound of the probability
of observing a clique of certain size and density in a random graph, which
appeared in (1):

$\displaystyle\mathcal{T}(\hat{B}_{c})=\exp\left(-\frac{1}{4}\zeta\left(\gamma_{1},p_{1}\right)|S_{c}||F_{c}|-\frac{1}{4}\zeta\left(\gamma_{2},p_{2}\right)|V_{c}|^{2}\right),$
(8)

where $\zeta(a,b)=\left\\{\frac{1}{(a-b)^{2}}+\frac{1}{3(a-b)}\right\\}^{-1}$,
$\gamma_{1}=\frac{\left|H_{c}\right|}{\left|S_{c}\right|\left|F_{c}\right|}$,
$p_{1}=\frac{\sum_{k\in S,(ij)\in F}I(\hat{\beta}_{(ij),k}\neq 0)}{|S||F|}$,
$\gamma_{2}=\frac{\left|F_{c}\right|}{\binom{|V_{c}|}{2}}$,
$p_{2}=\frac{\sum_{i,j\in V,i<j}I(\hat{\beta}_{(ij),k}\neq
0)}{\binom{|V|}{2}}$. We formally present our proposed network-based
permutation test for the significance of each extracted $\hat{B}_{c}$ in
Algorithm 2. The permutation procedure outlined in Algorithm 2 is effective in
simulating the null distribution of the test statistic $T(\hat{B}_{c})$.
Therefore, FWER can be controlled effectively, yielding a corrected $p$-value
for each extracted $\hat{B}_{c}$.

Algorithm 2 Assess the significance of $\\{\hat{B}_{c}\\}_{c=1,\ldots,C}$

Input: $\\{a_{(ij),k}\\},\hat{C}\geq 1,\\{\hat{B}_{c}\\},\alpha$

output: FWER-controlled significance values for each $\hat{B}_{c}$

1\. Choose a sound cut-off $\hat{r}$ and binarize graph
$B[\hat{r}]$:$B[\hat{r}]_{(ij),k}=I(a_{(ij),k}>\hat{r})$

2\. Estimate edge densities for $B$ and $G$, respectively:
$\displaystyle\text{Overall density}:\hat{p}_{1}=\frac{\sum_{k\in S,(ij)\in
F}I(a_{(ij),k}>\hat{r})}{|S||F|},\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{p}_{2}=\frac{\sum_{i,j\in V,i<j}f_{ij}}{\binom{|V|}{2}};$
$\displaystyle\text{Within-subgraph
density}:\hat{\gamma}_{1}=\frac{\sum_{k\in\hat{S}_{c},(ij)\in\hat{F}_{c}}I(a_{(ij),k}>\hat{r})}{|\hat{S}_{c}||\hat{F}_{c}|},\leavevmode\nobreak\
\leavevmode\nobreak\
\hat{\gamma}_{2}=\frac{\sum_{i,j\in\hat{V}_{c},i<j}f_{ij}}{\binom{|\hat{V}_{c}|}{2}}.$

3\. Calculate test statistic $\mathcal{T}_{0}(\hat{B}_{c})$ for the current
$\hat{B}_{c}$

4\. Shuffle group labels of data $P$ times and implement MOAT on each shuffled
graph

5\. Store the maximal test statistic for each simulation $l=1,\ldots,L:$
$\displaystyle\mathcal{T}_{l}=\sup_{c=1,\ldots,\hat{C}^{l}}\exp\left(-\frac{1}{4}\zeta\left(\gamma_{1}^{l},p_{1}^{l}\right)|\hat{S}^{l}_{c}||\hat{F}^{l}_{c}|-\frac{1}{4}\zeta\left(\gamma_{2}^{l},p_{2}^{l}\right)|\hat{V}^{l}_{c}|^{2}\right)$

6\. Calculate the percentile of $\mathcal{T}_{0}(\hat{B}_{c})$ in
$\\{\mathcal{T}_{l}\\}_{l\in(1,\ldots,L)}$ as the FWER $q$-value and reject
$\mathbb{H}_{0}$ if $q<\alpha$, where $\alpha$ is a pre-specified significant
level

Evaluating the joint effect of multiple SIs on FCs. With each
$\hat{B}_{c}=(S_{c},F_{c})$, we have a set of structural measures
$\\{X_{k},k\in S_{c}\\}$ associated with functional measures
$\\{Y_{ij},(ij))\in F_{c}\\}$. However, this does not automatically provide
the joint effect (i.e., $\sum_{(ij))\in F_{c},k\in S_{c}}\beta_{(ij),k}X_{k}$)
of the selected SIs on each selected FC measure. To assess the joint effect,
we can adopt the existing multivariate-to-multivarite analysis tools (e.g.,
CCA). Detailed procedures for applying CCA on $\hat{B}_{c}$ is provided in
Appendix D. Alternatively, one can conduct low-rank regression on outcomes and
predictors in each $\hat{B}_{c}$ to estimate the final effect size (Vounou et
al., , 2010; Wang et al., , 2012; Zhu et al., , 2014; Kong et al., , 2019).

##  3 Simulation

In this simulation study, we probed whether MOAT can extract informative
subgraphs $\\{\hat{B}_{c}\\}$ from the multi-level graph $B$ and $G$ with high
accuracy and replicability. We evaluate MOAT to finite-sample simulation data
under various conditions (e.g., different sample sizes and effect sizes) with
comparisons to several commonly used biclustering methods and sCCA-based
methods.

###  3.1 Synthetic data

We generate synthetic FC data
$\mathbf{Y}^{(d)}=\Big{\\{}y^{(d)}_{ij}\Big{\\}}_{i,j<n}$ and SI data
$\boldsymbol{X}^{(d)}=\big{(}x_{1}^{(d)},\ldots,x_{m}^{(d)}\big{)}^{T}$ based
on the following multivariate Gaussian distribution:

$\displaystyle\begin{pmatrix}\boldsymbol{X}^{(d)}\\\
\mathbf{Y}^{(d)}\end{pmatrix}$ $\displaystyle\sim
N\begin{bmatrix}\begin{pmatrix}\boldsymbol{\mu}_{X}\\\
\boldsymbol{\mu}_{Y}\end{pmatrix}\\!\\!,&\left(\begin{array}[]{cc}\mathbf{\Sigma}_{X,X}&\mathbf{\Sigma}_{X,Y}\\\
\mathbf{\Sigma}_{Y,X}&\mathbf{\Sigma}_{Y,Y}\end{array}\right)\end{bmatrix},$
(9)

where $\begin{pmatrix}\boldsymbol{\mu}_{X}\\\
\boldsymbol{\mu}_{Y}\end{pmatrix}$ is the partitioned mean vector of SI and FC
data respectively, and
$\mathbf{\Sigma}=\left(\begin{array}[]{cc}\mathbf{\Sigma}_{X,X}&\mathbf{\Sigma}_{X,Y}\\\
\mathbf{\Sigma}_{Y,X}&\mathbf{\Sigma}_{Y,Y}\end{array}\right)$ is the
partitioned variance-covariance matrix. For simplicity, we set
$\begin{pmatrix}\boldsymbol{\mu}_{X}\\\ \boldsymbol{\mu}_{Y}\end{pmatrix}$ as
a zero vector, representing normalized data, while the construction of
$\mathbf{\Sigma}$ depends on two key factors: the multi-level network
structure and effect sizes (i.e., FC-SI association strength). Both of these
factors are elaborated upon in the following paragraph.

In this simulation, we consider $500$ SIs and $4950$ FCs, where the FC
measures are calculated based on a brain network with $100$ regions, resulting
in ${100\choose 2}=4950$ pairwise connectivity values. To determine the
network structure for $\mathbf{\Sigma}$, we consider the following multi-level
graph consisting of (i) a bipartite graph $B=(S,F;H)$ depicting the FC-SI
associations, where $|S|=500,|F|=4950$; (ii) a unipartite graph $G=(V;E)$
depicting the brain functional connectome, where $|V|=100$. Specially, we
generate two sub-networks within $B$, denoted as $B_{1}$ and $B_{2}$,
characterized by higher FC-SI partial correlations $\rho_{1}$ and $\rho_{2}$
than the rest of $B$. $B_{1}$ consists of $40$ SI measures and $435$ FC
measures, where the $435$ FC measures collectively compose a functional
connectome $G_{1}$ of $30$ brain regions; $B_{2}$ consists of $60$ SI measures
and $190$ FC measures, where the $190$ FC measures collectively compose
another functional connectome $G_{2}$ of $20$ brain regions. For a visual
representation of these two sub-networks, please refer to the graph
illustration in Figure 3.

Figure 3: Applying MOAT and comparative methods to an inference matrix
$\boldsymbol{A}_{500\times 4950}$ with sample size $D=200$.
$\boldsymbol{A}_{500\times 4950}$ is simulated to store partial correlations
between $500$ SI measures and $4950$ FC measures, where the $4950$ FC measures
are obtained from region-wise functional connections between 100 brain
regions. Two sub-networks $B_{1}$ and $B_{2}$ with FC-SI partial correlations
$(\rho_{1},\rho_{2})=(0.40,0.35)$ are generated within $B$, while the partial
correlations for the rest of $B$ are set to be $\rho_{0}=0.15$. The right
panel shows the results of applying MOAT, biclustering and sCCA methods. Both
MOAT and biclustering methods are more accurate in revealing sub-network
patterns than sCCA methods since they incorporate cluster/network information.
When jointly evaluating TPR and TNR, MOAT outperforms other methods, as TPR
and TNR rely heavily on accurate sub-network extraction and inference.

Built on this network architecture, we configure the covariance matrix
$\mathbf{\Sigma}$ such that $\rho_{1},\rho_{2}>\rho_{0}$ to emulate different
effect sizes. Here, we set $\rho_{0}=0.15$ as the partial correlation of FC-SI
edges outside of $B_{1}$ and $B_{2}$. Next, by correlating the FC and SI data
simulated from (9) using the aforespecified
$\begin{pmatrix}\boldsymbol{\mu}_{X}\\\ \boldsymbol{\mu}_{Y}\end{pmatrix}$ and
$\mathbf{\Sigma}$, we obtain an FC-SI association matrix
$\boldsymbol{A}_{500\times 4950}$. $\boldsymbol{A}$ governs the edge variable
$H$ in the bipartite graph $B=(S,F;H)$ by $h_{(ij),k}=I(a_{(ij),k}>r)$, where
$r$ is a pre-selected threshold for correlation strength. Lastly, to assess
MOAT performance under different settings, three configurations of
$(\rho_{0},\rho_{1},\rho_{2};D)$ are simulated: $(0.15,0.55,0.60;200)$,
$(0.15,0.60,0.45;300)$, and $(0.15,0.70,0.40;400)$, where $D$ represents the
sample size as defined previously. For each configuration, we simulate 500
repeated data sets $\\{\boldsymbol{A}^{l}\\}_{l\in(1,\ldots,500)}$ to better
access accuracy and replicability of MOAT.

Figure 4: Table1:  Inference results of MOAT and comparative methods under
different settings. $\rho_{1},\rho_{2}$ denote the partial correlations
between SI and FC measures in the dense sub-networks $B_{1}$ and $B_{2}$,
respectively, and $D$ denotes the sample size. We summarize the means
(standard deviations) of the TPR and TNR for identified SIs, FCs, and FC-SI
pairs based on 500 repeated simulations. The results show that MOAT
outperforms the compared methods under different scenarios, indicating its
effectiveness in identifying the correct dense sub-networks with true positive
FC-SI associated pairs.

###  3.2 Performance evaluation

For each simulated dataset, we apply MOAT to estimate the multi-level sub-
networks $\hat{B}_{c}$ containing strong FC-SI associations and perform our
proposed network-based permutation test outlined in Algorithm 2 on
$\boldsymbol{A}^{l}$. Regarding ${B}_{c}$ extraction and
$\\{\beta_{(ij),k}\neq 0\\}$ identification, we benchmark MOAT against a few
popular appoaches including (i) three biclustering methods that are commonly
used for sub-network detection: Bipartite Spectral Graph Partitioning (BSGP)
(Wieling and Nerbonne, , 2009), Information Theoretic Learning (ITL)
(Erdogmus, , 2002), and Factor Analysis for Bicluster Information Acquisition
(FABIA) (Hochreiter et al., , 2010); (ii) two sCCA-based methods that identify
and measure the associations between two canonical/latent types of variables:
a Large-Scale Sparse Kernel Canonical Correlation method proposed by Uurtio et
al.,  (2019), and sCCA through a penalized matrix decomposition (sCCA-PMA)
proposed by Witten et al.,  (2009).

We evaluate methods’ performance by assessing the deviation of the estimated
$\hat{B}_{c}$ from true $B_{c}$ at both node-level, and edge-level (i.e.,
$\hat{\beta}_{(ij),k}\neq 0$ v.s true $\beta_{(ij),k}\neq 0$). Specifically,
we consider the comparions from the following three perspectives: SI variable
selection, FC variable selection, and FC-SI pair selection. We use true
positive rate (TPR) and true negative rate (TNR) as the evaluation criteria
for both node-level and edge-level deviations. TPR is determined by the
proportion of FC/SI nodes or FC-SI edges in $B_{c}$ that can be recovered by
$\hat{B}_{c}$; TNR is determined by the proportion of FC/SI nodes or FC-SI
edges in $B/B_{c}$ that can be recovered by $B/\hat{B}_{c}$.

Figure 3 provides a graphical overview of the performance of each method.
Table 1 demonstrates the performance of all methods under multiple settings.
The TPR and TNR are determined by the accuracy of both sub-network extraction
and network-level inference. In general, both MOAT and biclustering-based
methods can recover sub-network patterns more accurately than sCCA-based
methods because the network structures of FC-SI association patterns can be
better recognized. Under different settings, MOAT can detect the target sub-
networks with high sensitivity with few or none false-positive FC-SI edges
because the cost of removing a true positive association or including a false
positive edge is very high, as regulated by objective function (5). The
performance of biclustering methods is also improved with increased effect
sizes with low false positive rates and medium to low sensitivity. In
contrast, sCCA-based methods is invariant to different effect sizes, and may
miss the underlying FC-SI association patterns due to various noise.

Overall, MOAT is robust to noise and sensitive to organized FC-SI association
patterns. MOAT outperforms comparable biclustering and sCCA methods under
different settings, especially when systematic FC-SI association patterns are
present. This superiority stems from MOAT’s ability to accurately extract FC-
SI association patterns through multi-level sub-network analysis and tailored
sub-network-level inference.

##  4 Study of FC-SI associations in brain connectome data

###  4.1 UK Biobank sample and neuroimaging data

We aim to investigate the systematic effects of certain structural brain
imaging measures on the functional connectome using UK Biobank data (Sudlow et
al., , 2015). The UK Biobank is a vast biomedical database with approximately
half a million participants from the UK, where a total of 40,923 healthy
individuals were found to have usable resting-state fMRI (rs-fMRI) data that
passed quality control (Alfaro-Almagro et al., , 2018). Among them, a subgroup
of 4,242 individuals possessed complete data on the following three sets of
measurements we have chosen to focus on in this study:

  1. (i)

105 SI measures: we collected 105 SI variables including 39 white matter
integrity measures and 66 cortical thickness measures. The white matter
integrity reflects the overall health and coherence of brain white matter and
was assessed by fractional anisotropy (FA) obtained from DTI data in this
study. The DTI data was pre-processed using ENIGMA DTI protocols (Jahanshad et
al., , 2013) and white matter tracts were labeled based on the JHU ICBM DTI-81
Atlas (Smith et al., , 2006; Mori et al., , 2008). A complete list of the 39
regional white matter tracts can be found in Appendix C.3. On the other hand,
cortical thickness measures gauge the width of the gray matter of the human
cortex, and were obtained from T1 MRI and labeled based on the FreeSurfer
atlas (Tustison et al., , 2014).

  2. (ii)

30,135 FC measures: functional connectome data were obtained from rs-fMRI data
based on Brainnetome Atlas (Fan et al., , 2016). We first performed rs-fMRI
preprocessing for all participants and then extracted the averaged time series
of blood-oxygen-level-dependent (BOLD) signals from 246 functional brain
regions, resulting in $\binom{246}{2}=30,135$ region-pair FC measures. Details
of imaging acquisition and fMRI preprocessing are provided in Appendix C.1.

  3. (iii)

4 confounding variables: we adjusted four confounding variables including age
(years: $61.46\pm 7.40$), sex (M/F: 2003/2239), educational level (years:
$17.37\pm 3.92$), and body mass index (BMI) ($kg/m^{2}:26.35\pm 4.30$). These
variables have been used in previous neuroimaging literature on studying brain
functional connectivity (Miller et al., , 2016; Alfaro-Almagro et al., , 2021;
Bischof and Park, , 2015; Agustí et al., , 2018).

###  4.2 Results

We applied MOAT to the multimodal imaging data from the qualified 4,242 UK
Biobank participants. First, we obtained the FC-SI association inference
matrix $\boldsymbol{A}_{105\times 30135}$. Each entry in $\boldsymbol{A}$ is
$a_{(ij),k}=-\log(p_{(ij),k})$, where $p_{(ij),k}$ represents the $p$-value
testing the association between the $k$-th SI measure and the FC outcome
between two brain regions $i$ and $j$. Next, we performed a hard-thresholding
sparsity constraint by setting $a_{(ij),k}=a_{(ij),k}I(a_{(ij),k}<\epsilon)$
for some positive integer $\epsilon$ (Zhang et al., , 2023). We then applied
our proposed greedy peeling algorithm 1 to the inference matrix
$\boldsymbol{A}$, with tuning parameters $\lambda_{1}=1.25,\lambda_{2}=1.5$
selected by the KL divergence with a mixed Bernoulli distribution based on
random graphs $B$ and $G$. Algorithm 1 returned one multi-level sub-network
$\hat{B_{1}}\in B$. Lastly, we performed the network-level statistical
inference on $\hat{B_{1}}$ using Algorithm 2. The testing results showed that
the systematic association pattern of $\hat{B_{1}}$ is statistically
significant ($p<0.0001$).

Figure 5: Application of MOAT on a real neuroimaging dataset obtained from the
UK Biobank. (a) shows a heatmap of the inference matrix
$\boldsymbol{A}_{105\times 30135}$. Each element is a $-log(p)$ value testing
the association between SI and FC measures. (b) shows a heatmap of a
reorganized $\boldsymbol{A}$ with element in the detected subgraph
$\hat{B}_{1}$ pushed to the top left corner. Within the second level $G$
composed by the selected FCs, there exists an organized subgraph $G_{1}\subset
G$ revealed in (e). As shown in (c) and (d), the FC-SI pairs within the
identified $\hat{B}_{1}$ have significantly stronger associations compared to
those outside of the network.

Specifically, results show that $\hat{B_{1}}$ comprised $|S_{1}|=23$ SI
measures and $|F_{1}|=1316$ FC outcomes, as highlighted in Figure 5(b).
Furthermore, the extracted $|F_{1}|$ unfolded into a dense clique
$\hat{G}_{1}\in G$ consisting of $|V_{1}|=79$ regions, as illustrated in
Figure 5(e). The FC-SI pairs within the identified sub-network $\hat{B}_{1}$
demonstrate significantly stronger associations compared to those outside of
the network, as evidenced by the high $R^{2}$ and $t$-statistics shown in
Figure 5 (c-d). The 23 extracted SI measures consist of 3 cortical thickness
measures and 20 FAs: the three cortical thickness measures correspond to the
mean thickness of the parahippocampal, superior temporal, and cuneus gyrus;
while for the 20 FA measures extracted, the top four with the strongest FC
associations are CST-R (corticospinal tract, right hemisphere), CST-L
(corticospinal tract, left hemisphere), ICP (inferior cerebellar peduncles),
and FX (fornix). More detailed information about the remaining 16 FA measures
can be found in Appendix E.3. Figure 6 (left panel) illustrates the names and
spatial locations of the 20 selected FAs.

The right panel in Figure 6 shows the spatial distributions of
within-$\hat{B}_{1}$ brain regions (79 regions in total), where they are
predominantly located in six cortices: frontal, subcortical, temporal,
parietal, insular, and limbic. Moreover, these regions consist of several
well-defined brain functional networks including temporo-frontal, somatomotor,
ventral attention, frontoparietal, and (partial) default mode network (DMN).
Overall, Figure 6 provides a 3D demonstration showcasing the systematic
association patterns between the subsets of SIs and FCs revealed by MOAT.
Notably, both the FC-SI significantly associated sub-network ($\hat{B}_{1}$)
and the brain functional sub-connectome ($\hat{G}_{1}$) exhibit well-organized
topological structures.

Figure 6: Extracted FC-SI associated sub-networks by MOAT, adjusted for age,
sex, educational level, and BMI. Specifically, the left panel shows the SI
cluster in the extracted sub-networks, which contains 20 FA measures and 3
cortical thickness measures. The top four FA measures with remarkably stronger
associations with the detected FC sub-network are CST-R, CST-L, ICP, and FX.
The three selected cortical thickness measures correspond to the mean
thickness of the parahippocampal, superior temporal, and cuneus gyrus. The
right panel shows the detected FC sub-network that is strongly impacted by the
23 selected SIs. The detected FC sub-network covers 79 brain regions and
composes several high-level cognitive brain networks including the default
mode network, temporo-frontal network, somatomotor, ventral attention, and
frontoparietal network.

Figure 7: Table 2:  20 selected white matter tracts strongly associated with
identified FC sub-network

We further applied CCA on the extracted sub-network $\hat{B}_{1}$ to
quantitatively measure the canonical associations among the FC-SI pairs within
$\hat{B}_{1}$. Results showed that the sample canonical correlations of the
first three canonical variate pair in $\hat{B}_{1}$ were $0.81$, $0.69$, and
$0.68$ respectively. In contrast, we performed sparse CCA proposed by Witten
et al.,  (2009) on the full graph $B$ and $G$, given the ultra high
dimensionality of data. This yielded sample canonical correlations of $0.18$,
$0.15$, and $0.14$ for the first three canonical variate pairs, respectively.
Notably, MOAT can better recognize the underlying large-scale FC-SI
association patterns and then provide an improved estimation of the
multivariate-to-multivariate association.

In summary, the application of MOAT helps to unfold the complex yet
systematical and strong interplay between subsets of structural and functional
measures of the human brain. Our findings suggest i) FC-SI associations are
highly concentrated in a subset of SIs and FC sub-networks rather than
exhibiting a whole-brain diffuse distribution pattern; ii) several FC sub-
networks are primarily influenced by white matter integrity measures (refer to
Table 2 in Appendix E.3 for possible mapping relationships); iii) multiple SI
measures jointly affect the overall FC outcomes based on MOAT-guided CCA
analysis. While, on a high level, our results align well with previous medical
findings (Cheung et al., , 2008; Pradat and Dib, , 2009; Chaddock-Heyman et
al., , 2013; Corrigan et al., , 2015), MOAT reveals more refined patterns with
improved spatial specificity and biological interpretability.

##  5 Discussion

Our newly developed approach, MOAT, offers a novel strategy to investigate the
complex association patterns between multimodal neuroimaging data with matrix-
outcomes (FCs) and a vector of imaging predictors (SIs). MOAT deciphers the
complex FC-SI association patterns in a multi-level graph structure revealing
the joint effect of a small set of SI predictors on FC sub-networks. The
multi-level graph structure can effectively reduce the number of parameters
while preserving the spatial specificity of FCs and SIs. MOAT delivers
findings in organized multi-level sub-networks largely suppressing individual
false positive FC-SI associations (see Lemma 1 in section 2.4). We developed
computationally efficient algorithms to extract multi-level sub-networks. We
further showed the consistency of the MOAT method. In addition, we develop a
tailored network-level inference approach to test the extracted multi-level
sub-networks while controlling FWER. Last, MOAT is also compatible with
existing multivariate-to-multivariate analysis tools (e.g., CCA).

In our case study, we investigated the FC-SI associations based on a large
sample and revealed systematic association patterns with neurological
explanations. This may enhance our understanding of how the brain structure
and function interactively work during resting states and may lead to insights
that can guide future cognitive and psychiatric therapy. However, since UK
biobank participants mainly consist of elder Caucasians, our conclusion may be
limited. Further investigation and integrated analysis is required to gain
more comprehensive understanding of the FC-SI associations. The software
package for MOAT is available at https://github.com/TongLu-
bit/MultilayerNetworks-MOAT.

Declaration of interest: none.

Acknowledgments

Tong Lu and Shuo Chen were supported by the National Institutes of Health
under Award Numbers 1DP1DA04896801, EB008432, and EB008281. Yuan Zhang was
supported by the National Science Foundation under Award Number DMS-2311109.

## References

  * Agustí et al., (2018) Agustí, A., García-Pardo, M. P., López-Almela, I., Campillo, I., Maes, M., Romaní-Pérez, M., and Sanz, Y. (2018).  Interplay between the gut-brain axis, obesity and cognitive function.  Frontiers in neuroscience, 12:155.
  * Alfaro-Almagro et al., (2018) Alfaro-Almagro, F., Jenkinson, M., Bangerter, N. K., Andersson, J. L., Griffanti, L., Douaud, G., Sotiropoulos, S. N., Jbabdi, S., Hernandez-Fernandez, M., Vallee, E., et al. (2018).  Image processing and quality control for the first 10,000 brain imaging datasets from uk biobank.  Neuroimage, 166:400–424.
  * Alfaro-Almagro et al., (2021) Alfaro-Almagro, F., McCarthy, P., Afyouni, S., Andersson, J. L., Bastiani, M., Miller, K. L., Nichols, T. E., and Smith, S. M. (2021).  Confound modelling in uk biobank brain imaging.  NeuroImage, 224:117002.
  * Bahrami et al., (2019) Bahrami, M., Laurienti, P. J., and Simpson, S. L. (2019).  Analysis of brain subnetworks within the context of their whole-brain networks.  Human brain mapping, 40(17):5123–5141.
  * Bai et al., (2009) Bai, F., Zhang, Z., Watson, D. R., Yu, H., Shi, Y., Yuan, Y., Qian, Y., and Jia, J. (2009).  Abnormal integrity of association fiber tracts in amnestic mild cognitive impairment.  Journal of the neurological sciences, 278(1-2):102–106.
  * Ball et al., (2017) Ball, G., Aljabar, P., Nongena, P., Kennea, N., Gonzalez-Cinca, N., Falconer, S., Chew, A. T., Harper, N., Wurie, J., Rutherford, M. A., et al. (2017).  Multimodal image analysis of clinical influences on preterm brain development.  Annals of neurology, 82(2):233–246.
  * Bischof and Park, (2015) Bischof, G. N. and Park, D. C. (2015).  Obesity and aging: Consequences for cognition, brain structure and brain function.  Psychosomatic medicine, 77(6):697.
  * Bowman et al., (2012) Bowman, F. D., Zhang, L., Derado, G., and Chen, S. (2012).  Determining functional connectivity using fmri data with diffusion-based anatomical weighting.  NeuroImage, 62(3):1769–1779.
  * Buckner et al., (2008) Buckner, R. L., Andrews-Hanna, J. R., and Schacter, D. L. (2008).  The brain’s default network: anatomy, function, and relevance to disease.  Annals of the new York Academy of Sciences, 1124(1):1–38.
  * Bullmore and Sporns, (2009) Bullmore, E. and Sporns, O. (2009).  Complex brain networks: graph theoretical analysis of structural and functional systems.  Nature reviews neuroscience, 10(3):186–198.
  * Cao et al., (2014) Cao, M., Wang, J.-H., Dai, Z.-J., Cao, X.-Y., Jiang, L.-L., Fan, F.-M., Song, X.-W., Xia, M.-R., Shu, N., Dong, Q., et al. (2014).  Topological organization of the human brain functional connectome across the lifespan.  Developmental cognitive neuroscience, 7:76–93.
  * Chachlakis et al., (2019) Chachlakis, D. G., Prater-Bennette, A., and Markopoulos, P. P. (2019).  L1-norm tucker tensor decomposition.  IEEE Access, 7:178454–178465.
  * Chaddock-Heyman et al., (2013) Chaddock-Heyman, L., Erickson, K. I., Voss, M. W., Powers, J. P., Knecht, A. M., Pontifex, M. B., Drollette, E. S., Moore, R. D., Raine, L. B., Scudder, M. R., et al. (2013).  White matter microstructure is associated with cognitive control in children.  Biological psychology, 94(1):109–115.
  * Chekuri et al., (2022) Chekuri, C., Quanrud, K., and Torres, M. R. (2022).  Densest subgraph: Supermodularity, iterative peeling, and flow.  In Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 1531–1555. SIAM.
  * Chen et al., (2023) Chen, S., Zhang, Y., Wu, Q., Bi, C., Kochunov, P., and Hong, L. E. (2023).  Identifying covariate-related subnetworks for whole-brain connectome analysis.  Biostatistics, page kxad007.
  * Cheung et al., (2008) Cheung, V., Cheung, C., McAlonan, G., Deng, Y., Wong, J., Yip, L., Tai, K., Khong, P., Sham, P., and Chua, S. (2008).  A diffusion tensor imaging study of structural dysconnectivity in never-medicated, first-episode schizophrenia.  Psychological medicine, 38(6):877–885.
  * Corrigan et al., (2015) Corrigan, F., Grand, D., and Raju, R. (2015).  Brainspotting: Sustained attention, spinothalamic tracts, thalamocortical processing, and the healing of adaptive orientation truncated by traumatic experience.  Medical Hypotheses, 84(4):384–394.
  * Craddock et al., (2013) Craddock, R. C., Jbabdi, S., Yan, C.-G., Vogelstein, J. T., Castellanos, F. X., Di Martino, A., Kelly, C., Heberlein, K., Colcombe, S., and Milham, M. P. (2013).  Imaging human connectomes at the macroscale.  Nature methods, 10(6):524–539.
  * Drevets et al., (2008) Drevets, W. C., Price, J. L., and Furey, M. L. (2008).  Brain structural and functional abnormalities in mood disorders: implications for neurocircuitry models of depression.  Brain structure and function, 213:93–118.
  * Erdogmus, (2002) Erdogmus, D. (2002).  Information theoretic learning: Renyi’s entropy and its applications to adaptive system training.  University of Florida.
  * Fan et al., (2016) Fan, L., Li, H., Zhuo, J., Zhang, Y., Wang, J., Chen, L., Yang, Z., Chu, C., Xie, S., Laird, A. R., et al. (2016).  The human brainnetome atlas: a new brain atlas based on connectional architecture.  Cerebral cortex, 26(8):3508–3526.
  * Fornito et al., (2016) Fornito, A., Zalesky, A., and Bullmore, E. (2016).  Fundamentals of brain network analysis.  Academic press.
  * Goeman et al., (2022) Goeman, J. J., Górecki, P., Monajemi, R., Chen, X., Nichols, T. E., and Weeda, W. (2022).  Cluster extent inference revisited: quantification and localization of brain activity.  arXiv preprint arXiv:2208.04780.
  * Hayden et al., (2006) Hayden, E. P., Wiegand, R. E., Meyer, E. T., Bauer, L. O., O’connor, S. J., Nurnberger Jr, J. I., Chorlian, D. B., Porjesz, B., and Begleiter, H. (2006).  Patterns of regional brain activity in alcohol-dependent subjects.  Alcoholism: Clinical and Experimental Research, 30(12):1986–1991.
  * Hochreiter et al., (2010) Hochreiter, S., Bodenhofer, U., Heusel, M., Mayr, A., Mitterecker, A., Kasim, A., Khamiakova, T., Van Sanden, S., Lin, D., Talloen, W., et al. (2010).  Fabia: factor analysis for bicluster acquisition.  Bioinformatics, 26(12):1520–1527.
  * Honey et al., (2010) Honey, C. J., Thivierge, J.-P., and Sporns, O. (2010).  Can structure predict function in the human brain?  Neuroimage, 52(3):766–776.
  * Hotelling, (1933) Hotelling, H. (1933).  Analysis of a complex of statistical variables into principal components.  Journal of educational psychology, 24(6):417.
  * Jahanshad et al., (2013) Jahanshad, N., Kochunov, P. V., Sprooten, E., Mandl, R. C., Nichols, T. E., Almasy, L., Blangero, J., Brouwer, R. M., Curran, J. E., de Zubicaray, G. I., et al. (2013).  Multi-site genetic analysis of diffusion images and voxelwise heritability analysis: A pilot project of the enigma–dti working group.  Neuroimage, 81:455–469.
  * Johnson and Sinanovic, (2001) Johnson, D. and Sinanovic, S. (2001).  Symmetrizing the kullback-leibler distance.  IEEE Transactions on Information Theory.
  * Jolliffe and Cadima, (2016) Jolliffe, I. T. and Cadima, J. (2016).  Principal component analysis: a review and recent developments.  Philosophical transactions of the royal society A: Mathematical, Physical and Engineering Sciences, 374(2065):20150202.
  * Kemmer et al., (2018) Kemmer, P. B., Wang, Y., Bowman, F. D., Mayberg, H., and Guo, Y. (2018).  Evaluating the strength of structural connectivity underlying brain functional networks.  Brain Connectivity, 8(10):579–594.
  * Kong et al., (2019) Kong, D., An, B., Zhang, J., and Zhu, H. (2019).  L2rm: Low-rank linear regression models for high-dimensional matrix responses.  Journal of the American Statistical Association.
  * Lasky-Su et al., (2008) Lasky-Su, J., Neale, B. M., Franke, B., Anney, R. J., Zhou, K., Maller, J. B., Vasquez, A. A., Chen, W., Asherson, P., Buitelaar, J., et al. (2008).  Genome-wide association scan of quantitative traits for attention deficit hyperactivity disorder identifies novel associations and confirms candidate gene associations.  American Journal of Medical Genetics Part B: Neuropsychiatric Genetics, 147(8):1345–1354.
  * Li et al., (2012) Li, Y., Long, J., He, L., Lu, H., Gu, Z., and Sun, P. (2012).  A sparse representation-based algorithm for pattern localization in brain imaging data analysis.  PloS one, 7(12):e50332.
  * Lin et al., (2013) Lin, D., Zhang, J., Li, J., Calhoun, V. D., Deng, H.-W., and Wang, Y.-P. (2013).  Group sparse canonical correlation analysis for genomic data integration.  BMC bioinformatics, 14(1):1–16.
  * Lu et al., (2023) Lu, T., Zhang, Y., Kochunov, P., Hong, E., and Chen, S. (2023).  Network method for voxel-pair-level brain connectivity analysis under spatial-contiguity constraints.  arXiv preprint arXiv:2305.01596.
  * Marek et al., (2022) Marek, S., Tervo-Clemmens, B., Calabro, F. J., Montez, D. F., Kay, B. P., Hatoum, A. S., Donohue, M. R., Foran, W., Miller, R. L., Hendrickson, T. J., et al. (2022).  Reproducible brain-wide association studies require thousands of individuals.  Nature, 603(7902):654–660.
  * Margulis and Sagan, (2000) Margulis, L. and Sagan, D. (2000).  What is life? Univ of California Press.
  * Mbatchou et al., (2021) Mbatchou, J., Barnard, L., Backman, J., Marcketta, A., Kosmicki, J. A., Ziyatdinov, A., Benner, C., O’Dushlaine, C., Barber, M., Boutkov, B., et al. (2021).  Computationally efficient whole-genome regression for quantitative and binary traits.  Nature genetics, 53(7):1097–1103.
  * Mihalik et al., (2022) Mihalik, A., Chapman, J., Adams, R. A., Winter, N. R., Ferreira, F. S., Shawe-Taylor, J., Mourão-Miranda, J., Initiative, A. D. N., et al. (2022).  Canonical correlation analysis and partial least squares for identifying brain-behaviour associations: a tutorial and a comparative study.  Biological Psychiatry: Cognitive Neuroscience and Neuroimaging.
  * Miller et al., (2016) Miller, K. L., Alfaro-Almagro, F., Bangerter, N. K., Thomas, D. L., Yacoub, E., Xu, J., Bartsch, A. J., Jbabdi, S., Sotiropoulos, S. N., Andersson, J. L., et al. (2016).  Multimodal population brain imaging in the uk biobank prospective epidemiological study.  Nature neuroscience, 19(11):1523–1536.
  * Mori et al., (2008) Mori, S., Oishi, K., Jiang, H., Jiang, L., Li, X., Akhter, K., Hua, K., Faria, A. V., Mahmood, A., Woods, R., et al. (2008).  Stereotaxic white matter atlas based on diffusion tensor imaging in an icbm template.  Neuroimage, 40(2):570–582.
  * Nichols, (2012) Nichols, T. E. (2012).  Multiple testing corrections, nonparametric methods, and random field theory.  Neuroimage, 62(2):811–815.
  * Pradat and Dib, (2009) Pradat, P.-F. and Dib, M. (2009).  Biomarkers in amyotrophic lateral sclerosis: facts and future horizons.  Molecular diagnosis & therapy, 13:115–125.
  * Schaid et al., (2018) Schaid, D. J., Chen, W., and Larson, N. B. (2018).  From genome-wide associations to candidate causal variants by statistical fine-mapping.  Nature Reviews Genetics, 19(8):491–504.
  * Smith et al., (2006) Smith, S. M., Jenkinson, M., Johansen-Berg, H., Rueckert, D., Nichols, T. E., Mackay, C. E., Watkins, K. E., Ciccarelli, O., Cader, M. Z., Matthews, P. M., et al. (2006).  Tract-based spatial statistics: voxelwise analysis of multi-subject diffusion data.  Neuroimage, 31(4):1487–1505.
  * Smith et al., (2004) Smith, S. M., Jenkinson, M., Woolrich, M. W., Beckmann, C. F., Behrens, T. E., Johansen-Berg, H., Bannister, P. R., De Luca, M., Drobnjak, I., Flitney, D. E., et al. (2004).  Advances in functional and structural mr image analysis and implementation as fsl.  Neuroimage, 23:S208–S219.
  * Sudlow et al., (2015) Sudlow, C., Gallacher, J., Allen, N., Beral, V., Burton, P., Danesh, J., Downey, P., Elliott, P., Green, J., Landray, M., et al. (2015).  Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  PLoS medicine, 12(3):e1001779.
  * Sun et al., (2022) Sun, D., Rakesh, G., Haswell, C. C., Logue, M., Baird, C. L., O’Leary, E. N., Cotton, A. S., Xie, H., Tamburrino, M., Chen, T., et al. (2022).  A comparison of methods to harmonize cortical thickness measurements across scanners and sites.  Neuroimage, 261:119509.
  * Tang et al., (2016) Tang, Y., Liu, X., Wang, J., Li, M., Wang, Q., Tian, F., Su, Z., Pan, Y., Liu, D., Lipka, A. E., et al. (2016).  Gapit version 2: an enhanced integrated tool for genomic association and prediction.  The plant genome, 9(2):plantgenome2015–11.
  * Tustison et al., (2014) Tustison, N. J., Cook, P. A., Klein, A., Song, G., Das, S. R., Duda, J. T., Kandel, B. M., van Strien, N., Stone, J. R., Gee, J. C., et al. (2014).  Large-scale evaluation of ants and freesurfer cortical thickness measurements.  Neuroimage, 99:166–179.
  * Uurtio et al., (2019) Uurtio, V., Bhadra, S., and Rousu, J. (2019).  Large-scale sparse kernel canonical correlation analysis.  In International Conference on Machine Learning, pages 6383–6391. PMLR.
  * Vounou et al., (2010) Vounou, M., Nichols, T. E., Montana, G., Initiative, A. D. N., et al. (2010).  Discovering genetic associations with high-dimensional neuroimaging phenotypes: A sparse reduced-rank regression approach.  Neuroimage, 53(3):1147–1159.
  * Wang et al., (2012) Wang, H., Nie, F., Huang, H., Kim, S., Nho, K., Risacher, S. L., Saykin, A. J., Shen, L., and Initiative, A. D. N. (2012).  Identifying quantitative trait loci via group-sparse multitask regression and feature selection: an imaging genetics study of the adni cohort.  Bioinformatics, 28(2):229–237.
  * Wang et al., (2011) Wang, H., Nie, F., Huang, H., Risacher, S., Ding, C., Saykin, A. J., and Shen, L. (2011).  Sparse multi-task regression and feature selection to identify brain imaging predictors for memory performance.  In 2011 International Conference on Computer Vision, pages 557–562. IEEE.
  * Wang et al., (2023) Wang, Y., Yan, G., Wang, X., Li, S., Peng, L., Tudorascu, D. L., and Zhang, T. (2023).  A variational bayesian approach to identifying whole-brain directed networks with fmri data.  The Annals of Applied Statistics, 17(1):518–538.
  * Wang et al., (2020) Wang, Z., Novikov, A., Zolna, K., Merel, J. S., Springenberg, J. T., Reed, S. E., Shahriari, B., Siegel, N., Gulcehre, C., Heess, N., et al. (2020).  Critic regularized regression.  Advances in Neural Information Processing Systems, 33:7768–7778.
  * Wehrle et al., (2020) Wehrle, F. M., Lustenberger, C., Buchmann, A., Latal, B., Hagmann, C. F., O’Gorman, R. L., and Huber, R. (2020).  Multimodal assessment shows misalignment of structural and functional thalamocortical connectivity in children and adolescents born very preterm.  Neuroimage, 215:116779.
  * Wieling and Nerbonne, (2009) Wieling, M. and Nerbonne, J. (2009).  Bipartite spectral graph partitioning to co-cluster varieties and sound correspondences in dialectology.  In Proceedings of the 2009 Workshop on Graph-based Methods for Natural Language Processing (TextGraphs-4), pages 14–22.
  * Wig et al., (2014) Wig, G. S., Laumann, T. O., and Petersen, S. E. (2014).  An approach for parcellating human cortical areas using resting-state correlations.  Neuroimage, 93:276–291.
  * Witten et al., (2009) Witten, D. M., Tibshirani, R., and Hastie, T. (2009).  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  Biostatistics, 10(3):515–534.
  * Woo et al., (2014) Woo, C.-W., Krishnan, A., and Wager, T. D. (2014).  Cluster-extent based thresholding in fmri analyses: pitfalls and recommendations.  Neuroimage, 91:412–419.
  * Wu et al., (2021) Wu, Q., Zhang, Y., Huang, X., Ma, T., Hong, L. E., Kochunov, P., and Chen, S. (2021).  A multivariate to multivariate approach for voxel-wise genome-wide association analysis.  bioRxiv, pages 2021–11.
  * Yohai, (2008) Yohai, V. J. (2008).  Optimal robust estimates using the kullback–leibler divergence.  Statistics & probability letters, 78(13):1811–1816.
  * Zalesky et al., (2010) Zalesky, A., Fornito, A., and Bullmore, E. T. (2010).  Network-based statistic: identifying differences in brain networks.  Neuroimage, 53(4):1197–1207.
  * Zhang and Xia, (2018) Zhang, A. and Xia, D. (2018).  Tensor svd: Statistical and computational limits.  IEEE Transactions on Information Theory, 64(11):7311–7338.
  * Zhang et al., (2023) Zhang, J., Sun, W. W., and Li, L. (2023).  Generalized connectivity matrix response regression with applications in brain connectivity studies.  Journal of Computational and Graphical Statistics, 32(1):252–262.
  * Zhang et al., (2022) Zhang, J., Wang, H., Zhao, Y., Guo, L., and Du, L. (2022).  Identification of multimodal brain imaging association via a parameter decomposition based sparse multi-view canonical correlation analysis method.  BMC bioinformatics, 23(3):1–14.
  * Zhao et al., (2023) Zhao, Z., Chen, C., Adhikari, B. M., Hong, L. E., Kochunov, P., and Chen, S. (2023).  Mediation analysis for high-dimensional mediators and outcomes with an application to multimodal imaging data.  Computational Statistics & Data Analysis, 185:107765.
  * Zhou and Li, (2014) Zhou, H. and Li, L. (2014).  Regularized matrix regression.  Journal of the Royal Statistical Society. Series B, Statistical Methodology, 76(2):463.
  * Zhu et al., (2014) Zhu, H., Khondker, Z., Lu, Z., and Ibrahim, J. G. (2014).  Bayesian generalized low rank regression models for neuroimaging phenotypes and genetic markers.  Journal of the American Statistical Association, 109(507):977–990.
  * Zhu et al., (2017) Zhu, X., Suk, H.-I., Wang, L., Lee, S.-W., Shen, D., Initiative, A. D. N., et al. (2017).  A novel relational regularization feature selection method for joint regression and classification in ad diagnosis.  Medical image analysis, 38:205–214.
  * Zhuang et al., (2017) Zhuang, X., Yang, Z., Curran, T., Byrd, R., Nandy, R., and Cordes, D. (2017).  A family of locally constrained cca models for detecting activation patterns in fmri.  NeuroImage, 149:63–84.