Functional connectivity across the human subcortical auditory system using an autoregressive matrix-Gaussian copula graphical model approach with partial correlations

The mammalian auditory pathway conveys acoustic information between the inner ear—where sounds are mechanoelectrically transduced in the cochlea—and auditory cortex by way of multiple subcortical nuclei across the brainstem, midbrain, and thalamus. Much of our knowledge of the auditory system arises from anatomical and physiological research with nonhuman animal models (Brugge, 1992; McIntosh & Gonzalez-Lima, 1991; Webster, 1992). This work has contributed tremendously to our understanding of the mammalian auditory system. However, due to methodological and ethical limitations, our ability to directly assess auditory function in the human nervous system is severely constrained (Moerel et al., 2021).

This is particularly true for subcortical auditory structures. In mammals, the ascending central auditory pathway receives signals from the cochlea of the inner ear by way of the cochlear nerve, which principally innervates the cochlear nucleus in the brainstem. Auditory signals are then transmitted to the superior olivary complex, which is the first decussation point at which signals largely pass contralaterally from the left cochlear nucleus to the right superior olive (and similarly from right to left). From the superior olive, auditory signals travel (via the lateral lemniscus) to the inferior colliculus in the midbrain. The last subcortical auditory structure is the medial geniculate nucleus of the thalamus, which then passes information to the primary auditory cortex. In addition to the ascending “lemniscal” auditory pathway, an equal number of efferent connections transmit top-down information from higher order auditory regions to earlier auditory structures (Malmierca & Ryugo, 2011; Winer, 2005). Due to the small size of the subcortical auditory structures—tightly packed in with other heterogeneous nuclei and white matter pathways—and their anatomical location deep within the cranium, the subcortical auditory structures have received limited attention in noninvasive human research.

Functional magnetic resonance imaging (fMRI) is the most popular noninvasive method for probing macroscopic network-related brain activity. While studies of the human subcortical auditory system are somewhat limited, previous task-based fMRI research has functionally localized the subcortical auditory structures (Sitek et al., 2019), identified the tonotopic frequency mappings within the auditory midbrain and thalamus (De Martino et al., 2013; Moerel et al., 2015; Ress & Chandrasekaran, 2013), separated top-down and bottom-up speech-selective subregions of auditory thalamus (Mihai et al., 2019; Tabas et al., 2021), and recorded level-dependent BOLD signals throughout the auditory pathway (Hawley et al., 2005; Sigalovsky & Melcher, 2006).

In contrast to task-related BOLD activity, fMRI connectivity methods (often utilizing resting-state fMRI paradigms and full correlation analysis) are commonly used to assay cortical brain networks (Biswal et al., 1995; Gordon, Laumann, Adeyemo, et al., 2017; Power et al., 2011; Smith, Beckmann, et al., 2013), including the cortical auditory system (Abrams et al., 2020; Cha et al., 2016; Chen et al., 2020; Eckert et al., 2008; Maudoux et al., 2012; Ren et al., 2021). However, fMRI connectivity methods have limited history in subcortical research, especially in the auditory system, where they have—to our knowledge—only been utilized a handful of times to assess connectivity differences between individuals with and without tinnitus percepts (Berlot et al., 2020; Hofmeier et al., 2018; Leaver et al., 2016; J. Zhang et al., 2015).

From the seminal resting-state connectivity studies identifying human default mode and motor networks (Biswal et al., 1995; Fox & Raichle, 2007), to work linking functional connectivity with behavioral variability (Baldassarre et al., 2012; Deng et al., 2016), to investigations into brain network differences associated with disorders (Chai et al., 2016; Di Martino et al., 2011; Ferri et al., 2018; Greicius et al., 2007; Hahn et al., 2011; Husain & Schmidt, 2014; Kaiser et al., 2015; Sitek et al., 2016; Wilson et al., 2022), full correlation analysis has contributed tremendously to our understanding of human brain networks. However, moving beyond the traditional full correlation analysis should enable greater specificity in assessing functional connectivity patterns, particularly in identifying specific node-to-node connectivity patterns within an established brain network (Marrelec et al., 2006; Smith, 2012). In contrast to full correlations, which represent both direct and indirect connections, partial correlation analyses represent the direct association between two specific nodes after filtering out the effects of the remaining nodes and thus hold great promise for estimating direct functional connectivity within a network (Smith, Beckmann, et al., 2013). (Refer to the illustration in Section S1.1 in the Supplementary Materials.) For these reasons, partial correlation approaches are increasingly used to study functional connectivity networks in the brain (Wang et al., 2016; Warnick et al., 2018), including improved identification of individualized connectivity profiles compared with full correlation methods (Menon & Krishnamurthy, 2019), as well as improved prediction of brain disorders (de Vos et al., 2018; Reeves et al., 2023; Skåtun et al., 2017) and identification of general cognitive ability (Sripada et al., 2021). However, we are unaware of the application of such methods to assess functional connectivity within subcortical networks, particularly within the human auditory system.

In this article, we build on the Bayesian precision factor model (PFM) introduced recently in Chandra et al. (2021) to develop a novel highly robust autoregressive matrix-Gaussian copula graphical model (ARMGCGM) to assess partial correlation-based functional connectivity in a specific network in the human brain that spans subcortical and cortical regions, the auditory system. Notably, partial correlations can be readily obtained from the precision matrix, i.e., the inverse of the covariance matrix. The PFM decomposes the model precision matrix into a flexible low-rank and diagonal structure, then exploits that to design efficient estimation algorithms. However, PFM makes the assumption that data are independent and identically distributed (iid) multivariate Gaussian random variables. Several studies in the literature also make this assumption (Yu et al., 2022; L. Zhang et al., 2014) which has the very useful implication that a zero partial correlation between two variables (equivalent to a zero entry in the precision matrix in the corresponding position) also means independence between them after removing the effects of other variables. However, in many applications—including ours—the iid Gaussian assumption can be restrictive as data are often not Gaussian distributed. Additionally, data from successive fMRI volumes exhibit strong autocorrelation. The ARMGCGM extends the PFM to the case where the univariate marginals can be any arbitrary distribution while also accounting for the autocorrelations between successive fMRI scans. The association between the variables is modeled using a Gaussian copula that implies conditional independence for zero partial correlation, allowing easy interpretability of the conditional dependence graph. As developed, the ARMGCGM approach is broadly applicable for studying undirected functional graphs using large-scale fMRI data.

We use the novel ARMGCGM to investigate functional connectivity between specific nodes across the human auditory system. We used publicly available 7T resting-state fMRI from over 100 participants to examine auditory connectivity. To probe connectivity within the auditory system, we included auditory cortical regions of interest as well as subcortical auditory regions derived from human histology (Sitek et al., 2019). As the auditory pathway comprises a chain of multiple subcortical structures, and due to the largely lateralized organization of the lemniscal auditory system, we hypothesized that partial correlations would be greatest between adjacent nodes in the same hemisphere, over and above the contributions from other auditory (and nonauditory) regions of interest. In particular, because of the critical position of auditory midbrain and thalamus as computational hubs involving bottom-up and top-down information transfer (Mihai et al., 2019; Tabas et al., 2021), we hypothesized strong partial correlations between inferior colliculus and medial geniculate. We further assessed reliability across acquisitions by separately analyzing data with anterior–posterior and posterior–anterior phase-encoding directions, as well as leave-10%-out cross-validation. We then compared our ARMGCGM-based connectivity results with those from a full correlation approach as well as alternative partial correlation methods. Overall, our consistent findings of hierarchical connectivity within the auditory system using our novel partial correlation method—consistent across data partitions and leave-10%-out validation— demonstrate the methodological reliability of our ARMGCGM approach as well as the neurobiological organization of auditory structures in the human primary auditory system.

We used resting-state fMRI from 106 participants in the 7T Human Connectome Project (Elam et al., 2021; Van Essen et al., 2012). Specifically, we utilized the minimally preprocessed volumetric data in common space (Glasser et al., 2013). BOLD fMRI data were acquired with 1.6 mm isotropic voxel size across four runs of a resting-state paradigm (repetition time [TR] = 1 s, 900 TRs per run). Two runs were acquired with anterior–posterior (AP) phase encoding, and two were acquired with posterior–anterior (PA) phase encoding. For all regions of interest (ROIs) and runs, we discarded the first 50 TRs to increase stability.

For each individual and each run, we extracted mean timeseries from predefined regions of interest (ROIs). Subcortical auditory ROIs were defined using the Sitek–Gulban atlas (Sitek et al., 2019). Cortical ROIs were defined using FreeSurfer’s implementation of the DKT atlas (Dale et al., 1999; Klein & Tourville, 2012). For this study we used transverse temporal gyrus (TTG) and superior temporal gyrus (STG) as auditory ROIs, as well as pericalcarine cortex (Calc) and superior frontal gyrus (SFG) as nonauditory control ROIs (see Section S1.5 and Fig. S3 in Supplementary Materials). Mean timeseries were extracted for each ROI using Nilearn’s [Masker] function. For an overview of materials and methods used in this study, please see Fig. S4 in Supplementary Materials Section 1.6.

BOLD fMRI is prone to geometric distortions in the phase-encoding direction which can be largely corrected using a variety of methods (Esteban et al., 2021; Jezzard & Balaban, 1995). Adjacent to motion-sensitive cerebrospinal fluid (CSF), the brainstem is particularly susceptible to such geometric distortions. Although the HCP minimal preprocessing pipeline corrects for phase-encoding distortions (Glasser et al., 2013), to isolate the potential residual contribution of phase-encoding direction on functional connectivity estimates, we conducted separate analyses on data collected with the posterior–anterior (PA) phase-encoding direction (runs 1 and 3) and the anterior–posterior (AP) phase-encoding direction (runs 2 and 4) and compared the results. As the fMRI data acquired in the two phases will be analyzed separately but using the same probability-model, we use the same notations for the different phases to describe our proposed model in the following sections.

We further evaluated our results using a leave-10%-out cross-validation (CV) approach. We assessed the stability of the connectivity estimates by checking the correlation between the results obtained from each acquisition scheme.

The Bayesian precision factor model (PFM), developed recently in Chandra et al. (2021), provides a novel computationally efficient robust technique for estimating precision matrices. Since partial correlations can be readily obtained from the precision matrix, the approach allows straightforward estimation of the underlying connectivity graphs. A plausible alternative could be to obtain the precision matrix by modeling the covariance matrix and then inverting its estimate (Lee & Kim, 2021). However, this approach often tends to exhibit poor empirical performance (Pourahmadi, 2013) and hence direct estimation of the precision matrix, if possible, may be preferable.

Most partial correlation-based conditional dependency graph estimation procedures in the statistical literature (Cai et al., 2020; Chandra et al., 2021; Friedman et al., 2008; Warnick et al., 2018) assume that the joint distribution of the data is the independent and identically distributed (iid) multivariate Gaussian, which implies that the univariate marginal distributions are also all Gaussians. However, successive fMRI scans have strong autocorrelations, and the marginal distributions for different ROIs exhibit substantial deviance from the Gaussian assumption. In this paper, we, therefore, extend the PFM to accommodate non-Gaussian marginals, while also appropriately accounting for their temporal dependence between successive scans.

Let ${Yt,j(r,i)}t=1T$ be the fMRI timeseries corresponding to the $i$-th individual’s $j$-th ROI at the $t$-th time point in the $r$-th run. In our application we are interested in studying the connectivity between $d=12$ ROIs along the central auditory pathway using fMRI timeseries of length $T=850$ from $N=106$ individuals each undergoing $R=2$ runs. We let $fj(r,i)$ be the (unknown) marginal density of $Yt,j(r,i)$ with corresponding cumulative distribution function (CDF) $Fj(r,i)$⁠.

Copulas provide a broadly applicable class of tools that allow the joint distribution of $Yt,j(r,i)$ to be flexibly characterized by first modeling the univariate marginals $fj(r,i)$ and then hierarchically modeling their dependencies by mapping the $Fj(r,i)(Yt,j(r,i))$’s to a joint probability space. For Gaussian copulas, this is done by setting $Zt,j(r,i)=Φ−1{Fj(r,i)(Yt,j(r,i))}$⁠, where $Φ(⋅)$ is the CDF of a standard Gaussian distribution. This implies marginally $Zt,j(r,i)∼N(0,1)$ for all $r,i,j,t$⁠, where $N(μ,σ2)$ denotes a univariate Gaussian distribution with mean $μ$ and variance $σ2$⁠. We let $Z(r,i)=((Zt,j(r,i)))T×d$ denote the matrix of fMRI signals corresponding to the $i$-th individual in the $r$-th run in the transformed Gaussian space. The Gaussian copula assumption on $Ft,j(r,i)(Yt,j(r,i))$’s implies that the joint distribution of $Z(r,i)$ is Gaussian as well. The dependencies between the $Zt,j(r,i)$’s are, therefore, characterized entirely by their correlations. Additionally, since the dependence relationships between the observed $Yt,j(r,i)$’s are modeled only through $Zt,j(r,i)$’s, these correlations also completely characterize the dependencies between the $Yt,j(r,i)$’s.

Note that probabilistic dependencies exist between $Zt,j(r,i)$’s across both $j$ and $t$⁠; dependence across $j$ incurs due to the interaction between the ROIs, whereas dependence across $t$ occurs due to the temporal dependence between successive fMRI scans. We let $RΩ$ denote the $d×d$ correlation matrix accounting for the dependence across the $d$ different ROIs across all runs and individuals. While our main interest lies in estimating these dependencies between the ROIs, it is also crucial to consider the temporal dependence in the $Yt,j(r,i)$’s. Let $Z.,j(r,i)=(Z1,j(r,i),…,ZT,j(r,i))T$ be the $j$-th column of $Z(r,i)$⁠, i.e., the timeseries corresponding to the $i$-th individual’s $j$-th ROI in the $r$-th run.

We develop our model in a hierarchical manner. To begin with, we use autoregressive (AR) processes of order $L$ to model higher-order temporal dependencies in the $Zt,j(r,i)$’s as

with $Zt−t′,j(r,i)=0$ if $t′≥t$⁠. We assume separate $(β,ς2)$ parameters across $(r,i,j)$ in equation (1) to make the model adapt to different timeseries patterns across different ROIs, individuals, and runs. Let $ξ.,j(r,i)$ be the standardized AR-corrected timeseries corresponding to the $j$-th ROI for $j=1,…,d$⁠. Define the $T×d$ matrix $Ξ(r,i)=(ξ⋅,1(r,i)…ξ⋅,d(r,i))=(ξ1,⋅(r,i)…ξT,⋅(r,i))T$ where $ξt,⋅(r,i)$ is the $t$-th row of $Ξ(r,i)$ and can be interpreted as the fMRI signals in the Gaussian copula space subsequent to filtering out the temporal dependence at time point $t$⁠. We then let

From the property of Gaussian copula, the probabilistic dependence structure between the ROIs is entirely characterized by the correlation matrix $RΩ$⁠. This strategy thus allows borrowing of information across subjects as well as runs in estimating the common correlation matrix $RΩ$ resulting in improved statistical precision. The formulations in equations (1)-(2) imply the following joint distribution on $Z(r,i)$

where $RT,j(r,i)$ is the correlation matrix of $Z⋅,j(r,i)$ implied by the $AR(L)$ model in equation (1).

Note that while the AR-corrected timeseries $ξ⋅,j(r,i)$’s are required to implement the model using MCMC, they can be obtained from the autoregressive parameters ${β1,j(r,i),…,βL,j(r,i),ςj2(r,i)}$ without having to invert the massive $(T−1)×(T−1)$ dimensional matrices $RT,j(r,i)$⁠. Exploiting the conditionally linear structure from equation (1), the $ξ⋅,j(r,i)$’s can be obtained without actually factorizing or inverting the $RT,j(r,i)$’s. Avoiding such expensive matrix inversion operations makes our model implementation highly scalable; we discuss the details in Section S1.2 of Supplementary Materials.

Next, let $Ω=((ωj,j′))$ be a precision matrix corresponding to $RΩ$⁠, i.e., $RΩ=Ψ−12Ω−1Ψ−12$⁠. with $Ψ=diag(Ω−1)$⁠. Then, by properties of copula and the simple multivariate Gaussian distribution, it can be shown that, irrespective of the form of the marginal $fj(r,i)$’s, the ROIs $Yt,j(r,i)$ and $Yt,j′(r,i)$ will be conditionally dependent on (and hence functionally connected with) each other given the rest if and only if $ωj,j′≠0$⁠. This way, the precision matrix $Ω$ characterizes the functional connectivity between the different $Yt,j(r,i)$’s across $j$⁠.

To estimate $Ω$⁠, we consider the framework of Chandra et al. (2021). We assume $Ω$ to admit a lower-rank plus diagonal (LRD) decomposition $Ω=ΛΛT+Δ$ where $Λ$ is a $d×q$ matrix and a diagonal matrix $Δ=diag(δ12,…,δd2)$ with positive $δj2$’s. Note that all positive definite matrices admit such a representation for some $q≤d$⁠.

Finally, we model the unknown univariate marginal distributions $fj(r,i)$’s using location-scale mixtures of Gaussians. Such mixtures can flexibly approximate a very large class of densities (Ghosal et al., 1999). Specifically, we let

where $πh,j(r,i)$ is the weight attached to the $h$-th mixture component and $∑h=1Kπh,j(r,i)=1$ with $K$ being a suitably chosen, moderately large, fixed integer.

While fMRI timeseries often showcase very distinct patterns and distributions for different individuals as well as across different ROIs of the same individual, we expect primary sensory region to be highly similar between healthy participants due to strongly stereotyped processing across individuals and well-conserved auditory function in these brain regions across evolution (Hutchison et al., 2013). In equation (4), we, therefore, consider separate mixture models across different ROIs, individuals, and runs, whereas in equation (3), the $Z(r,i)$’s share a common correlation matrix $RΩ$ across all individuals and runs. This modeling strategy also allows borrowing of information across individuals and runs to amplify signals for estimating resting-state functional connectivity in adult humans while accounting for subject and run-specific variabilities using separate marginal distributions. This is conceptually similar to approaches used in group independent component analysis (Calhoun et al., 2001) and cohort-level brain mapping (Varoquaux et al., 2013). In later sections, we showed that this approach yields highly consistent estimates of functional connectivity graphs even though BOLD signals from small deep brain regions have very low signal-to-noise ratio (Bianciardi et al., 2016; Colizoli et al., 2020; de Hollander et al., 2017; Sclocco et al., 2018).

With fMRI timeseries data ${Yt,j(r,i)}t=1T$ corresponding to the $i$-th individual’s $j$-th ROI at the $t$-th time point in the $r$-th run, our hierarchical approach can be summarized as the following:

• The marginal densities $fj(r,i)$ of $Y1,j(r,i),…,YT,j(r,i)$ are modeled using flexible location-scale mixtures of univariate Gaussians as specified in equation (4).

• Letting $Fj(r,i)$ be the corresponding CDF of $Y1,j(r,i),…,YT,j(r,i)$⁠, the timeseries are transformed into a Gaussian space $Zt,j(r,i)=Φ−1{Fj(r,i)(Yt,j(r,i))}$ implying that marginally $Zt,j(r,i)∼N(0,1)$ for all $r,i,j,t$⁠.

• Autocorrelation-corrected $ξt,j(r,i)$’s are obtained from the $Zt,j(r,i)$’s using the $AR(L)$ model from equation (1). Accordingly, the $ξt,.(r,i)=(ξt,1,…,ξt,d)T​′s$ are iid $N(0,RΩ)$ random vectors across all $r,i,t$⁠.

• Let $RΩ=Ψ−12Ω−1Ψ−12$ with $Ψ=diag(Ω−1)$⁠, where $Ω$ is the precision matrix of the $ξt,.(r,i)$’s. The precision factor model framework of Chandra et al. (2021) is then used to estimate $Ω$⁠.

In Figure 1 we provide a graphical illustration of the hierarchical structure of the ARMGCGM.

We take a Bayesian route to estimation and inference, where we assign priors to the model parameters, and then infer them based on samples drawn from the posterior using a Markov chain Monte Carlo (MCMC) algorithm discussed in detail in the Supplementary Materials. As was shown in Chandra et al. (2021), the LRD representation makes the MCMC sampling very efficient via a latent variable augmentation scheme. In the Supplementary Materials, we also discuss a multiple hypothesis testing-based edge discovery procedure that utilizes the posterior uncertainty of the parameters and controls the false discovery rate (FDR) at 10% level in a principled manner.

For the autoregressive parameters in equation (1), for all $r,i,j$ we assume

$β1,j(r,i)​,…,βL,j(r,i)​|ςj2(r,i)∼iidN(0,νβ−1ςj2(r,i)),ςj−2(r,i)∼iidGa(aς,bς),$where $νβ,aς,bς>0$ are fixed hyperparameters, and $Ga(a,b)$ denotes a gamma distribution with mean $a/b$ and variance $a/b2$⁠. For the parameters of the mixture models specifying the marginals in equation (4), we consider the following priors

where $(μ,σ2)∼NIG(μ0,ν0,a0,b0)$ is the normal-inverse-gamma prior implying that $μ|σ2∼N(μ0,ν0−1σ2)$ and $σ−2∼Ga(a0,b0)$⁠.

Ideally, $(β1,j(r,i),…,βL,j(r,i),ςj(r,i))$ should satisfy some constraints to conform with the $Zt,j(r,i)$’s being stationary and having unit variance. Enforcement of such constraints, however, presents a significant challenge. Instead, in our implementation, we left these parameters unrestricted but observed the posterior to always converge to regions where such properties were accurately satisfied.

Next, we consider a shrinkage prior on the elements of $Λ$ that shrinks redundant elements of $Λ$ to zero, allowing additional model-based parameter reduction. In particular, we assign a two-parameter generalization of the Dirichlet–Laplace (DL) prior from Bhattacharya et al. (2015) that allows more flexible tail behavior on the elements of $Λ$⁠. On a $d$-dimensional vector $θ=(θ1,…,θd)T$⁠, our DL prior with parameters $a$ and $b$⁠, denoted by $DL(a,b)$⁠, can be specified in the following hierarchical manner.

where $θj$ is the $j$-th element of $θ$⁠, $ϕ$ and $ϱ$ are vectors of same length as $θ$⁠, $Exp(a)$ is an exponential distribution with mean $1/a$⁠, and $Dir(a1,…,ad)$ is a $d$-dimensional Dirichlet distribution. The original DL prior is a special case with $b=1/2$⁠. We let $vec(Λ)∼DL(a,b)$⁠.

We use a Dirichlet process (DP) prior (Ferguson, 1973) on the $δj2$’s as $δj2|G∼iidG,G|α∼DP(α,G0) with G0=$$Ga(aδ,bδ),α∼Ga(aα,bα),$ where $α$ is the concentration parameter and $G0$ is the base measure to favor a smaller number of unique $δj2$’s in a model-based manner. The DP model allows clustering the $δj2$’s facilitating additional data-adaptive parameter reduction when necessary. Additionally, a DP prior has full prior support on the range space of the number of unique $δj2$’s implying a fully flexible model (see Chapter 4 of Ghosal and van der Vaart (2017)).

We discuss a Markov chain Monte Carlo (MCMC)-based strategy of sampling from the posterior of ARMGCGM in Section S1.2 of the Supplementary Materials, where the involving steps are parallelized over the subjects, allowing scalability. In Section S1.3 of the Supplementary Materials, we discuss the choice of hyperparameters used for the analyses presented in this paper. Our implementation using the proposed parallelized MCMC scheme ran in 125 min in a system with 13th Gen Intel(R) Core(TM) i9-13900K CPU and 128 GB RAM, fitting the ARMGCGM to the 12-node auditory network with 7,500 MCMC iterations, including both phase-encoding schemes.

Using our ARMGCGM approach, we first estimate the precision matrix $Ω$ and subsequently compute the partial correlation matrix. We report the significant edges subject to controlling the posterior false discovery rate (Sarkar et al., 2008) at the 10% level. (Details are provided in Section S1.4 in the Supplementary Materials.) In Figure 2 we provide the circos plots of the connectivity graphs along with respective weighted adjacency matrices. The $(j,h)$-th off-diagonal element of the adjacency matrices admit the value 0 if the edge between ROIs $j$ and $h$ is not statistically significant, if the edge is significant then the corresponding partial correlation is plugged in to indicate the strength of the edge.

We found that the strongest auditory connectivity was between adjacent structures in the same hemisphere, particularly between the auditory midbrain (inferior colliculus, or IC) and thalamus (medial geniculate body, or MGB) and between core and associative auditory cortex (TTG and STG). Minimal connectivity was observed between homologous auditory structures across hemispheres. Connectivity was also present between adjacent brainstem auditory structures (cochlear nucleus [CN] and superior olivary complex [SOC]), largely bilaterally. Interestingly, despite the SOC being the primary (and earliest) decussation point in the primary auditory pathway, we only observed partial correlations between right CN and left SOC (in both data partitions), not left CN and right SOC. (See Section 4 for potential explanations.)

Due to the anatomical location of the subcortical auditory structures—in dense, heterogeneous subcortical regions and largely adjacent to CSF—we conducted connectivity analyses separately on AP- and PA-acquired fMRI runs. We then compared the connectivity results from the two acquisition schemes. Overall connectivity patterns were highly similar between the two phase-encoding schemes, as seen in panels A and B of Figure 2. Subcortical connectivity was quite robust between the brainstem auditory regions. To quantify the similarity between the results in PA and AP acquisitions (plotted in Fig. 3A), we computed the Pearson correlation coefficient $r$ between the estimated partial correlations. We found a Pearson correlation $r=0.940$ with $p$-value $<0.001$ between the partial correlations for PA and AP acquisitions that were computed with the proposed ARMGCGM. Note that a large proportion of points overlap at the origin [0.0, 0.0], indicating that ARMGCGM induced sparsity in the connectivity graphs similarly in the two acquisition schemes, which contributes to the high Pearson correlation between the schemes. To assess the similarity in sign of connectivity between acquisition schemes, we computed the Jaccard dissimilarity on the signed off-diagonals of each adjacency matrix. For the proposed ARMGCGM, the Jaccard dissimilarity was 0.231. Additionally, we computed the Euclidean distance between partial correlations in the two acquisition schemes, compared it with some selected approaches from the literature and report the results in Table 1. These results indicate high level of consistency between the findings from the two different acquisition schemes, with the ARMGCGM exhibiting the strongest similarity.

To assess the stability of our proposed ARMGCGM, we ran leave-10%-out cross-validation by first randomly splitting the $N$ subjects into 10 equal folds. Then, for each of the 10 folds, we held out the BOLD signals ${Yt,j(r,i)}t=1T$ for all $i$’s in that fold and fit ARMGCGM on the remaining 90% subjects. Finally, we compared the estimated partial correlations across the folds. Similar approaches are used in the literature to check the stability of graph estimates (Gan et al., 2018; Mazumder & Hastie, 2012). We repeated this leave-10%-out cross-validation approach for both phase-encoding schemes. Partial correlation coefficients for each of the 20 folds (10 folds for each of the 2 acquisition schemes) are presented in Figure 3B, which shows highly similar partial correlation estimates across all folds and both phase-encoding schemes (intraclass correlation of edgewise partial correlations across cross-validation folds = 0.991).

In Section S1.7 of the Supplementary Materials, we performed additional simulation studies (Supplementary Fig. S5) which illustrates our method’s ability to efficiently recover the “true” underlying connectivity graph.

We compared with three standard approaches in the literature: (1) correlation analysis between the ROIs (Biswal et al., 1995; Cordes et al., 2000; Fox & Raichle, 2007; Lowe et al., 1998; Smith, Vidaurre, et al., 2013); and partial correlation analyses using (2) the graphical lasso (Friedman et al., 2008) and (3) the precision factor model (Chandra et al., 2021). In all comparisons, we did separate analyses for each of the acquisition schemes.

Here we study the marginal correlation between the ROIs. Letting $ρj,j′$ denote the correlation between ROIs $j$ and $j′$ in resting state we test

For each acquisition scheme, we first concatenate the timeseries across all runs and individuals. Then we perform $t$-tests for correlation for all $(j,j′)$ pairs followed by the Benjamini–Hochberg false discovery rate (FDR) correction (Benjamini & Hochberg, 1995) for multiplicity adjustment and control the FDR at level $0.10$⁠. Traditionally, full correlations are used to measure functional “connectivity” in resting-state fMRI studies (Biswal et al., 1995; Cordes et al., 2000; Lowe et al., 1998). In Figure 4, we provide the correlation graphs and correlation matrices separately for each acquisition scheme. We find the correlation graphs to be much denser compared with the partial correlation graphs presented in Figure 2. Unlike in the ARMGCGM method, we observed negative values when using full correlations (particularly in the PA acquisition scheme), although these negative correlations are generally close to 0. In this full correlation approach, the similarity between runs split by data acquisition scheme was characterized by Pearson’s $r$ of $0.811$ (⁠$p$-value $<0.001$⁠) and a Euclidean distance of 1.016. The Jaccard dissimilarity of the signed adjacency matrix was 0.338.

We first consider the graphical lasso (Glasso) approach (Friedman et al., 2008) as another alternative choice for partial correlation-based conditional graph estimation. Glasso assumes iid data from a multivariate Gaussian distribution (i.e., without any correction for autocorrelation) with l_1$(⋅)$ penalty on the precision matrix. In this analysis, we concatenated the $Yt,j(r,i)$ values across $(r,i)$ and created a $(TR)×d$ matrix, say $Y$⁠, for each acquisition scheme, and applied the Glasso model on $Y$⁠. We use 10-fold cross-validation to choose the optimal penalty parameter. The Glasso approach provides a point estimate of the sparse precision matrix and hence the functional connectivity network. In panels A and B of Figure 5, we plot the connectivity graphs and weighted adjacency matrices, respectively, separately for each acquisition scheme.

These functional connectivity graphs in Figure 5 are much denser compared with the estimates obtained by our proposed ARMGCGM in Figure 2. To quantify the robustness and stability of the Glasso graphs in this application, we compute Pearson correlation coefficient and Euclidean distance between the estimated partial correlations, and the Jaccard dissimilarity between the signs of the adjacency matrices in PA and AP acquisitions in the same manner as we did for ARMGCGM elaborated in Section 3.2. We reported the values in Table 1. Figure 5 indicates that the strong positive correlations are consistent across the acquisition schemes. However, substantial discrepancy can be observed for the weak edges, particularly for the negative partial correlations. This is a common phenomenon for graphical models if the Gaussian assumption is made on non-Gaussian distributed data; see, e.g., Section 5, example 1(a) in Guha et al. (2020).

In this analysis, we concatenated the $Yt,j(r,i)$ values across $(r,i)$ and created a $(TR)×d$ matrix, say $Y$⁠, for each acquisition scheme, and applied the PFM on $Y$⁠. Like Glasso, the PFM also assumes iid data from a multivariate Gaussian distribution and does not correct for temporal autocorrelations in the data. We infer on the graph using the Bayesian multiple comparison technique described in the Supplementary Materials. Results are provided in Figure 6. The connectivity graphs majorly differed with our ARMGCGM results, with the PFM-derived graph being much denser and includes more (weakly) negative edges. In the top panel of Figure 7, we plot the marginal Gaussian fits on the histograms of some BOLD signals. Clearly the simple Gaussian assumption does not hold here and a more sophisticated approach like ours is required.

To assess the robustness and stability of the PFM, we computed Pearson correlation coefficient and Euclidean distance between the estimated partial correlations, and the Jaccard dissimilarity between the signs of the adjacency matrices in PA and AP acquisitions in the same manner as we did for ARMGCGM elaborated in Section 3.2. We reported the values in Table 1, where we see that PFM exhibits more consistency than Glasso but ARMGCGM performed best.

We studied the goodness-of-fit of the proposed ARMGCGM. In the top panel of Figure 7, we plotted the sample histograms and the corresponding fitted marginal densities for some selected ROIs. The sample histograms strongly indicate the marginal distributions of the data to deviate substantially from Gaussian distributions, including some with multiple well-separated modes. Figure 7 also shows that our flexible location-scale mixture of Gaussians fits the data very well, even for the most complicated distributions.

Mixtures of Gaussians with reasonably large number of mixture components can approximate large classes of widely varying distributions. To validate whether the number of mixtures $(K=20)$ in model (4) is adequate, we computed the median number of nonempty clusters across MCMC samples for each subject and ROI in each run. In the bottom panel of Figure 7, we plotted the histograms of the medians across the subjects. As the number of nonempty clusters is smaller than $K$ consistently across all setups yielding excellent fits for complicated distributions, we conclude that our model specifications are adequate.

We recall from Section 2.3 that $Zt,j(r,i)$’s were the transformed BOLD signal timeseries in the Gaussian space and $ξt,j(r,i)$’s were the autocorrelation-corrected timeseries. To check for autocorrelation corrections using model (1), we plotted the partial autocorrelations between $Zt,j(r,i)$’s and $ξt,j(r,i)$’s across time for all ROIs. Since the $Z$ and $ξ$ values vary across MCMC iterations, we used their posterior expectations in this analysis. In Figure 8, we plot the partial autocorrelations for a randomly selected individual in run 1 of posterior-to-anterior acquisition scheme; the blue horizontal dashed lines represent the 5% band. The figure clearly indicates that the autoregressive model of order $L=5$ corrects for the autocorrelations. As we obtained very similar results for other runs and acquisition schemes, we have omitted them from the paper.

The mammalian auditory pathway consists of a series of obligatory and interconnected subcortical and cortical brain structures. Assessing the connectivity of the human subcortical auditory structures has been limited due to methodological challenges of noninvasive imaging of the deep, small structures. Recent acquisition and analytical advances enable finer grained investigations of connectivity throughout the brain, including the brainstem. In this paper, we validated a novel autoregressive matrix-Gaussian copula graphical model to estimate functional auditory connectivity patterns from a publicly available high-resolution resting-state functional MRI dataset. Using partial correlations (as opposed to full correlations) allowed us to identify specific relationships between nodes in a connectivity graph by removing shared variance across nodes (Supplementary Figs. S1 and S2). We found highly consistent connectivity patterns between adjacent auditory brain regions along the auditory pathway that demonstrate the efficacy of our connectivity method as well as the potential for functional connectivity investigations of the subcortical auditory system. Below, we separately discuss our novel scientific findings and our novel contributions to the statistics literature.

To date, there have been only limited applications of functional MRI methods to study subcortical auditory connectivity (Berlot et al., 2020; Hofmeier et al., 2018). Using our novel ARMGCGM approach in the present study, we found strong partial correlations between cochlear nucleus (CN) and superior olivary complex (SOC) bilaterally using resting-state functional MRI data. Most interestingly, we observed contralateral CN–SOC connectivity between right CN and left SOC (and in both data acquisition schemes), fitting the ground truth primary auditory pathway crossing from left to right (and vice versa) between CN and SOC (Barnes et al., 1943; Schofield, 1994). These functional connectivity patterns between CN and SOC have not been previously observed in human auditory brainstem in vivo but follow our understanding of the mammalian primary auditory pathway based on research in animal models (Barnes et al., 1943; Doucet & Ryugo, 2003; Harrison & Irving, 1966). Principally, auditory information that is transduced by the cochlea of each ear is transmitted via the cochlear nerve to the cochlear nucleus, the first stage of the central auditory pathway, on each side of the brainstem. In the primary auditory pathway, the lemniscal anteroventral subdivision of the cochlear nucleus enhances the fine temporal precision of incoming auditory signals (Pickles, 2015). From there, auditory signals are passed to both the ipsilateral and contralateral SOC for further auditory processing, including spatial localization (Moore, 2000). The SOC comprises multiple distinct subdivisions, which receive ipsilateral and contralateral connections from cochlear nucleus to varying degrees (Pickles, 2015), aligning with our overall bilateral connectivity results between CN and SOC.

While we observed consistent right CN and left SOC connectivity, it is unclear why similar patterns were not observed between left CN and right SOC. One contributor is the lower signal-to-noise ratio in fMRI data from the lower brainstem. Paired with the small size of each of the brainstem auditory nuclei, we may still be at the edge of what is detectable using present functional connectivity methods. Additionally, this analysis was conducted on “resting-state” fMRI data, during which no auditory stimuli of interest was presented or overt tasks were conducted. Resting-state fMRI connectivity in the cochlear nucleus and superior olivary complex has not been examined in the previous literature to our knowledge; it is possible that sound-evoked connectivity methods would evoke greater functional connectivity, particularly in these earliest stages of the auditory pathway. Further, ipsilateral connections (i.e., between left CN and left SOC and between right CN and right SOC) may be artifactually stronger due to their physical proximity. Even with relatively high 1.05 mm spatial resolution 7T fMRI data, CN and SOC on each side are only separated by a few voxels. These regions are thus at increased likelihood of sharing temporal fluctuations due to partial volume effects, wide point-spread functions, spatial dependence, or another as-yet-unsolved fMRI confounds that are particularly acute in the lower brainstem.

Moving up the primary auditory pathway, we observed significant partial correlations between ipsilateral inferior colliculus (IC) in midbrain and medial geniculate body (MGB) of the thalamus in both hemispheres and in both phase-encoding schemes. Inferior colliculus is a major convergence point in the auditory system, with the lemniscal IC subdivision being thought to convert distinct auditory features into discrete auditory objects for the first time in the auditory pathway. MGB continues the refinement of auditory objects via direct lemniscal connections from IC as well as rich corticofugal connections from auditory cortex to nonlemniscal MGB subdivisions (Pickles, 2015). Interestingly, we found strong partial correlations between left and right MGB in both datasets. Although not directly connected by large white matter bundles, left and right MGB are expected to process auditory information from IC at similar levels of abstraction. Thus, partial correlations may reflect indirect but shared neural mechanisms of auditory processing in the thalamus.

We did not observe IC partial correlation connectivity with either brainstem or cortical structures. IC is a key hub in the auditory system, receiving bottom-up sensory information as well as top-down modulating signals from auditory cortex and other brain regions. The lack of partial correlation connectivity with IC may be due to strong IC subdivision-specific functionality, with IC core primarily serving an ascending lemniscal role and dorsal and external IC having top-down and nonlemniscal functions. Averaging over these subdivisions may obfuscate specific connectivity patterns. Alternatively, our results may suggest that it does not have a specialized relationship with any one region beyond MGB but rather integrates and transforms auditory and other neural signals.

Finally in auditory cortex, transverse temporal gyrus (TTG)—the location of primary auditory cortex—was strongly connected with ipsilateral superior temporal gyrus (STG), which contains secondary and associative auditory cortices. Primary auditory cortex receives direct input from lemniscal MGB and is the last auditory structure with fine-grained tonotopicity (Pickles, 2015). In humans, STG is hierarchically structured, with portions further away from primary auditory cortex having increasingly wider temporal integration windows (Hamilton et al., 2018; Norman-Haignere et al., 2022) and greater categorical specificity (Bhaya-Grossman & Chang, 2022; Feng et al., 2021; Hamilton et al., 2020; Keshishian et al., 2023; Nourski et al., 2018; Pernet et al., 2015; Rauschecker & Tian, 2000; Rupp et al., 2022). While invasive recordings from human STG suggest a potential direct connection between MGB and posterior STG (Hamilton et al., 2021), we found mixed evidence for such a direct pathway in our partial correlation data (in one hemisphere in only one of the data partitions). Our partial correlation functional connectivity results align with a vast literature demonstrating information flow between primary and nonprimary auditory cortex (for review, see Hackett (2011). The lack of contralateral partial correlation approach measures the correlation between the concerned ROIs after filtering out the indirect effects of the remaining ROIs. Complementarily, some literature suggest that left and right auditory cortices process auditory information at distinct timescales and levels of abstraction (Güntürkün et al., 2020; Hickok & Poeppel, 2007; Zatorre et al., 2002), with the left auditory cortex being uniquely tuned to rapidly changing temporal information—such as the phonetics of speech sounds—while the right auditory cortex is more sensitive to slower changes in the spectral domain, particularly for speech prosody as well as music.

In comparison with our ARMGCGM partial correlation approach, we computed full correlations in the same network. Connections were much denser in the full correlation approach, aligning with the rich interconnectedness of the auditory system (Pickles, 2015). Unlike with partial correlations, which highlighted hierarchical connections between adjacent nodes along the auditory pathway, we observed positive full correlations between all auditory cortical regions, regardless of hemisphere. Additionally, we found a strong positive subnetwork including IC and MGB bilaterally, whereas many of these connections (such as between left and right IC) were absent in the ARMGCGM partial correlation analysis. Since partial correlations characterize connectivity between two nodes after filtering out the effects of the other nodes, our combined results point to widespread shared information across the auditory system (per full correlation analysis) with additional shared processing between adjacent nodes of the canonical auditory pathway (per partial correlation analysis). This suggests distinct but complementary use of full and partial correlations, with full correlation analysis identifying a rich network of interconnected nodes, while partial correlations are sensitive to strong node-to-node connections.

In this article, we developed an autoregressive matrix-Gaussian copula graphical model (ARMGCGM) for non-Gaussian distributed data with temporal autocorrelation, the problem of estimating brain connectivity patterns from resting-state fMRI data being the primary motivation. In the ARMGCGM, we begin by employing flexible individual-specific and brain-region-specific location-scale mixtures of Gaussians to model the marginal distributions of the observed fMRI data. Subsequently, these models are utilized to transform the observed data into latent Gaussian timeseries. We then apply higher-order autoregressive models to capture the temporal dependence in this transformed space within each region. Finally, a Gaussian copula is employed on the autocorrected timeseries to capture the conditional dependencies between different brain regions shared across individuals. The ARMGCGM thus allows borrowing of information across different subjects to infer on a common connectivity graph while taking into account subject and run-specific variability via flexible mixture models. We leverage recent advances on modeling precision matrices via a flexible but computationally efficient low-rank-diagonal decomposition method that not allows efficient exploration of the posterior space for estimating the connectivity graph. Compared with alternative approaches to exploiting partial correlations to estimate connectivity graphs, our proposed ARMGCGM method produces results that are more consistent across fMRI data acquisition schemes with respect to multiple metrics. Additionally, the results remain highly stable in a leave-10%-out cross-validation. Considering the low signal-to-noise ratio in BOLD signals from deep small auditory structures, our results demonstrate the sensitivity and specificity of our model to neurobiologically plausible connections. While the proposed ARMGCGM offers a more sophisticated method for modeling fMRI data, with the potential to fit arbitrarily complex marginal distributions, it comes with a higher numerical cost compared with simplistic parametric models. However, our parallelized Markov chain Monte Carlo implementation runs across all participants in just over 2 hours, showcasing a manageable computation time considering the dataset’s scale. In our view, the enhanced inference capabilities outweigh the computational burden.

Our study is the first to systematically assess connectivity across the human auditory pathway using multiple connectivity measures, with previous subcortical connectivity studies limited to full correlation analysis (Berlot et al., 2020; Hofmeier et al., 2018; Leaver et al., 2016; J. Zhang et al., 2015). Given the anatomical and methodological constraints with subcortical fMRI, the limited fMRI connectivity literature is not too surprising. First, the deep location of the subcortical structures places them far from MRI transmit and receive coils, limiting the signal-to-noise ratio from these regions (Miletić et al., 2020). Because the brainstem is relatively centrally located relative to the multiple receiver coils, accelerated acquisition techniques that are based on phase differences between receiver coils are less effective (Preibisch et al., 2015). Second, subcortical nuclei can be quite small, requiring higher resolution imaging protocols (which unfortunately trade off SNR in order to achieve greater spatial resolution). Third, subcortical nuclei are densely organized adjacent to nuclei with heterogeneous functions, so voxels immediately next to those containing core auditory structures could contain visual, motor, or sensory nuclei, white matter, CSF, or a combination of any of these. Ultimately, each of these constraints limits the SNR from subcortical auditory nuclei.

Constraints in human subcortical auditory research have translational consequences beyond basic science. For instance, while cochlear implants have been widely successful at providing sensory information to individuals with sensorineural hearing impairments with an intact cochlear nerve (Kral et al., 2019; Reiss, 2020), auditory prostheses in the central auditory system have been less successful (Lim et al., 2009; Shetty et al., 2021), due in no small part to our limited understanding of the complexity of sound representation in the ascending auditory pathway.

As the canonical neuroanatomy of the primary auditory pathway is consistent across individuals and well described in the literature (Sitek et al., 2019), we have the a priori expectation of a shared auditory graph across all participants. As the goal of this study was to map a network that is strongly expected based on anatomy and nonhuman neurophysiology, we built a joint model that includes data from all participants. This is similar to approaches used in group independent component analysis (Calhoun et al., 2001) and cohort-level brain mapping (Varoquaux et al., 2013). However, subject-specific differences in the distributions of the BOLD signals as well as autocorrelations between successive scans can induce artifactual and noisy edges in functional connectivity graphs, as seen in the Glasso and PFMs; we take care of these issues in ARMGCGM. Nevertheless, we did not investigate differences in functional connectivity between participants in the current article. Building on ARMGCGM to explore how functional connectivity varies between individuals and groups or as a function of behavior is a priority for future work.

In general, resting-state fMRI connectivity measures become more reliable with longer scans (Zuo et al., 2019). Measurement correlations increase as time in the scanner increases, from Pearson’s $r=0.82$ at 9 min to $r=0.92$ at 27 min to $r=0.97$ at 90 min (Laumann et al., 2015). Others described improved intraclass correlation coefficients with datasets beyond 20 minutes and up to 50 minutes (Xu et al., 2016). The Midnight Scan Club group (Gordon, Laumann, Gilmore, et al., 2017) computed a range of network connectivity metrics and found that reliability generally required at least 30 minutes of resting-state data per subject. One paper (Greene et al., 2020) specifically investigated functional connectivity in subcortical structures and found even longer scan requirements (up to 100 minutes) for subcortical structures due to decreased signal-to-noise ratios deeper in the brain. Additionally, primary sensory networks are among the most stable within and across participants (Gratton et al., 2018; Hutchison et al., 2013). We, therefore, believe it is appropriate and necessary to use datasets with longer scans of resting-state data in order to investigate even static subcortical auditory connectivity.

Further, while many brain networks exhibit temporal dynamics that can tell us about mental state (Fornito et al., 2012) or disease ((Sakoğlu et al., 2010); see Hutchison et al. (2013) for a review), functional connectivity within primary sensory networks is among the most stable over time (Gratton et al., 2018), as they share bidirectional physical connections, share contributions to the same physiological tasks, and are evolutionarily conserved across species (Hutchison & Everling, 2012). In the present work, we were interested in characterizing the stationary connectivity in the primary auditory pathway that is present at rest across individuals. Adapting time-varying dynamics into this model is a promising future direction, particularly if we are interested in higher level cognitive brain networks that vary as a function of task or mental state. Additionally, exploring recent advancements in the statistics literature (Ni et al., 2020; Song et al., 2020) to scale up MCMC implementations of Bayesian mixture models, exploring rank (Hoff, 2007), and empirical likelihood (Mengersen et al., 2013)-based methods for marginal estimation, etc. is of future interest. Finally, adapting our proposed undirected connectivity methods to estimate directed graphs holds promising neurobiological value and is a priority for ongoing and future research directions.

In this article, we validated a novel autoregressive matrix-Gaussian copula graphical model for partial correlation estimation while appropriately correcting for temporal autocorrelations. Using this approach, we identified functional connectivity in the human auditory system using resting-state functional MRI. Whereas a complementary approach using full correlations identified a rich network of interconnected auditory regions, partial correlations highlighted direct connections between adjacent structures along auditory pathways. In particular, subcortical connectivity was highly consistent across acquisitions, demonstrating the utility and applicability of functional connectivity methods in deep brain structures. In the future, we plan to investigate whole-brain partial correlation connectivity across sensory, motor, and higher cognitive networks using the proposed models and their relationship with behavior across individuals.

Data usage adhered to the Open Access Data Use Terms from the Human Connectome Project (WU-Minn HCP).

Data are publicly available through the Human Connectome Project (https://www.humanconnectome.org/study/hcp-young-adult). Analysis code is available in the Supplementary Materials.

Noirrit Kiran Chandra: Conceptualization (equal), Formal analysis (lead), Software (lead), Visualization (equal), Writing—original draft (equal), Writing—review and editing (equal). Kevin R. Sitek: Conceptualization (equal), Formal analysis (equal), Visualization (equal), Writing—original draft (lead), Writing—review and editing (equal). Bharath Chandrasekaran: Conceptualization (equal), Writing—Review and editing (equal). Abhra Sarkar: Conceptualization (equal), Writing—review and editing (equal).

This work was supported by the National Science Foundation grant NSF DMS-1953712 (to A.S. and B.C.) and the National Institutes of Health/National Institute on Deafness and Other Communication Disorders (K01DC019421 to K.R.S. and R01DC01550 to B.C.).

The authors have no competing interests.

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657), funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research, and by the McDonnell Center for Systems Neuroscience at Washington University.

Supplementary material for this article is available with the online version here: https://doi.org/10.1162/imag_a_00258.