Voxel-wise intermodal coupling analysis of two or more modalities using local covariance decomposition.

When individual subjects are imaged with multiple modalities, biological information is present not only within each modality, but also between modalities - that is, in how modalities covary at the voxel level. Previous studies have shown that local covariance structures between modalities, or intermodal coupling (IMCo), can be summarized for two modalities, and that two-modality IMCo reveals otherwise undiscovered patterns in neurodevelopment and certain diseases. However, previous IMCo methods are based on the slopes of local weighted linear regression lines, which are inherently asymmetric and limited to the two-modality setting. Here, we present a generalization of IMCo estimation which uses local covariance decompositions to define a symmetric, voxel-wise coupling coefficient that is valid for two or more modalities. We use this method to study coupling between cerebral blood flow, amplitude of low frequency fluctuations, and local connectivity in 803 subjects ages 8 through 22. We demonstrate that coupling is spatially heterogeneous, varies with respect to age and sex in neurodevelopment, and reveals patterns that are not present in individual modalities. As availability of multi-modal data continues to increase, principal-component-based IMCo (pIMCo) offers a powerful approach for summarizing relationships between multiple aspects of brain structure and function. An R package for estimating pIMCo is available at: https://github.com/hufengling/pIMCo.

INTRODUCTION

There is increased availability of multi‐modality neuroimaging data on individual subjects, with each modality containing unique information about brain structure or function. Such data allow us to explore patterns in individual modalities as well as patterns in the relationships between modalities, which we call intermodal coupling (IMCo), at global, regional, or local resolutions (Baller et al., 2022; Gu et al., 2021; Honey et al., 2009; Shokri‐Kojori et al., 2019; Tak et al., 2015; Uddin, 2013; Valcarcel, Linn, Khalid, et al., 2018a; Valcarcel, Linn, Vandekar, et al., 2018b; Vandekar et al., 2016). The progress made by these IMCo studies have transformed our understanding of the brain, and it suggests that advancements in the methodology for studying IMCo have the potential to further enable such insights.

On the global scale, intermodal relationships have long been of interest. For example, structural connectivity (SC) and functional connectivity (FC) are globally correlated in adults, but the relationship is less straightforward in children (Uddin, 2013). Regional studies have built on this global understanding. For example, studies have found SC‐FC coupling changes regionally in neurodevelopment to support cognition, and this coupling continues to evolve in aging, decreasing in sensorimotor systems while persevering in higher‐order cognitive systems (Baum et al., 2020; Esfahlani et al., 2021). Additionally, a number of studies have investigated how SC‐FC coupling varies across regions, showing this coupling to be strongest in primary sensory and motor regions, but that structure and function uncouple in higher‐level regions (Gu et al., 2021; Preti & Van De Ville, 2019; Vazquez‐Rodriguez et al., 2019).

In addition to structure–function relationships, regional relationships between metabolism and brain function have also been explored. In the study of energy utilization in the brain, Shokri‐Kojori et al. showed that regional correspondence between cerebral glucose metabolism and fluctuations in blood oxygenation not only differed between brain networks in healthy patients but was also sensitive to differences between patients with acute or chronic alcohol use (Shokri‐Kojori et al., 2019). Of note, these relationships were not identifiable by looking at individual modalities alone. Another regional study on metabolism‐function coupling identified significant associations between cerebral blood flow (CBF) and strength of functional connectivity in default, frontoparietal, and primary sensory‐motor networks. No significant association was found between CBF and functional connectivity strength in regions outside of these networks (Tak et al., 2015).

On the sub‐regional local scale, studies from our group on coupling between cortical thickness and sulcal depth have suggested the cortical sheet is generally thinner in sulcal locations than in gyral locations and that this relationship was more spatially heterogeneous than previously described (Vandekar et al., 2016). A separate study exploring local IMCo between CBF and resting‐state amplitude of low‐frequency fluctuations (ALFF) showed that age‐related declines in this measure of neurovascular coupling were most pronounced during mid‐adolescence and were enriched in the dorsal attention network (Baller et al., 2022). There were also differences in CBF‐ALFF coupling between males and females which were enriched in the frontoparietal network.

In these local IMCo studies, each vertex‐wise coupling value was defined as the slope of the weighted linear regression (WLR) best‐fit line for that local neighborhood between two modalities. However, because this method of calculating IMCo is based on regression slopes, it does not consider vertex‐level correlation and suffers from inherent asymmetry, where coupling values depend on which modality is defined as the independent variable in the WLR. This asymmetry necessitates arbitrary, yet influential, choices when it comes to analysis and limits straightforward interpretation. This measure for IMCo is also limited to only two modalities, so the study of coupling between more than two modalities using this method would require analysis of all pairwise couplings. As the number of total modalities increases, this approach can become challenging to interpret. Additionally, analysis of all pairwise couplings may not parsimoniously describe the overall degree of coupling across all modalities.

In response to these limitations, we propose a principal component analysis (PCA) based method for estimating IMCo that uses local covariance decomposition to define symmetric voxel‐wise coupling values valid for two or more modalities. This method reduces complex local covariance structures into a single value, thus providing an easily interpretable value that characterizes the strength of coupling in settings with two modalities. It also allows for simplified study of more complex local covariance structures in settings with more than two modalities.

To demonstrate its sensitivity to biologically relevant patterns, we show that PCA‐based IMCo (pIMCo) uncovers differences in three‐modality coupling between CBF, ALFF, and regional homogeneity (ReHo) with respect to age and sex in youth. We chose these modalities because local cortical coupling between vascular organization and spontaneous resting state fluctuations has been previously thought of as a measure of neurovascular coupling. Since ReHo describes local connectivity, or the degree of local neuronal activity coordination, we were interested in understanding the overall coupling between cerebral blood flow and both spontaneous and locally‐coordinated neuronal activity, which can be thought of as a more comprehensive measure of neurovascular coupling.

MATERIALS AND METHODS

Subjects

We included 803 subjects (340 males) from ages 8–22 (mean = 15.6; sd = 3.3) from the Philadelphia Neurodevelopmental Cohort (PNC) (Satterthwaite, Elliott, et al., 2014a). Of the 1601 PNC subjects who underwent neuroimaging, health screening as well as automated and manual image quality screening were performed. We excluded subjects in the following order: low T1‐weighted MRI quality (n = 61), low resting‐state fMRI (rfMRI) quality (n = 450), and low arterial spin labeling (ASL) quality (n = 54). Of the remaining subjects, we then excluded those meeting any of the following health exclusion criteria (n = 205): history of psychoactive medication, history of inpatient psychiatric hospitalization, or history of medical disorders that could impact brain function. Finally, ASL scans for which high‐quality partial volume correction could not be performed were excluded (n = 28). This resulted in the final set of 803 subjects used for this study.

The Institutional Review Boards of the University of Pennsylvania and the Children's Hospital of Pennsylvania approved all study procedures. All adult study subjects gave written informed consent; for subjects under the age of 18, parents or guardians provided written informed consent and subjects provided assent. Additional details of the PNC study have been previously described (Gur et al., 2020; Satterthwaite, Elliott, et al., 2014a).

Image acquisition

All PNC imaging was acquired using a single 3 T Siemens Tim Trio scanner with a 32‐channel head coil. To minimize motion, subjects' heads were stabilized using one foam pad over each ear and one foam pad over the top of the head. Image acquisition procedures have been previously described (Gur et al., 2020; Satterthwaite, Elliott, et al., 2014a).

T1‐weighted structural images were used for alignment of all scans into a common space. T1‐weighted images were acquired using a 3D‐encoded magnetization‐prepared, rapid‐acquisition gradient echo (MPRAGE) T1‐weighted sequence with the following settings: = 1810 ms; = 3.51 ms; FoV = 180 × 240 mm; matrix size = 192 x 256; number of slices = 160; slice thickness = 1 mm; inter‐slice gap = 0 mm; resolution = 0.9375 × 0.9375 × 1 mm. Cerebral blood flow (CBF) was estimated from a pseudo‐continuous arterial spin labeling (pcASL) sequence with a spin‐echo echoplanar readout and the following settings: = 4000 ms; = 15 ms; FoV = 220 × 220 mm; matrix size = 96 x 96; number of slices = 20; slice thickness = 5 mm; inter‐slice gap = 1 mm; resolution = 2.3 × 2.3 × 6 mm; 80 volumes. Maps of amplitude of low frequency fluctuations (ALFF) and regional homogeneity (ReHo) were estimated from 6 min of task‐free functional data from a blood‐oxygen‐level‐dependent (BOLD) weighted 2D EPI sequence with the following settings: = 3000 ms; = 32 ms; FoV = 192 × 192 mm; matrix size = 64 × 64; number of slices = 46; slice thickness = 3 mm; inter‐slice gap = 0 mm; resolution = 3 mm isotropic; 124 volumes. Subjects were instructed to stay awake, keep their eyes open, fixate on a displayed fixation cross, and remain still.

Image processing

Image processing of T1‐weighted structural images, pcASL scans, and rfMRI scans have been previously described (Baller et al., 2022; Gur et al., 2020). They are summarized here in brief. T1‐weighted structural images were processed using tools from Advanced Normalization Tools (ANTs) (Tustison et al., 2014). pcASL and rfMRI scans were processed using an eXtensible Connectivity Pipeline (XCP) which included tools from FSL and AFNI (Ciric et al., 2018; Cox, 1996; Jenkinson et al., 2012).

CBF was quantified from control‐label pairs using ASLtbx as previously described in (Satterthwaite et al. 2014b; Wang et al., 2008). Briefly, this quantification involved measuring the differences in signal between control and label acquisitions and then using the standard kinetic model to calculate the CBF estimate. Partial volume correction was performed on these CBF images using Bayesian Inference for Arterial Spin Labeling MRI (BASIL) (Chappell et al., 2009; Chappell et al., 2011).

For rfMRI processing, the XCP pipeline included: (1) field inhomogeneity correction with FSL FUGUE, (2) removal of initial rfMRI volumes, (3) alignment of volumes within the time series to a selected reference volume using FSL MCFLIRT, (4) interpolation of intensity outliers with AFNI 3dDespike, and (5) demeaning and removal of linear or quadratic trends. Images were then denoized using a 36‐parameter confound regression model that has been shown to minimize impact of motion artifact (Ciric et al., 2017). Finally, BOLD‐weighted time series as well as artifactual model time series were filtered using a first‐order Butterworth filter with a passband between 0.01 and 0.08 Hertz.

Voxel‐wise ALFF was defined as the sum of frequency bins between 0.01 and 0.08 Hertz using a Fourier transform of the time‐domain signal (Yang et al., 2007). Voxel‐wise ReHo was defined as Kendall's coefficient of concordance computed over the rfMRI time series in each voxel's 26‐voxel local neighborhood (Zang et al., 2004). Voxel‐wise maps were smoothed with a 6 mm full width at half maximum (FWHM) kernel to improve signal‐to‐noise ratio. CBF, ALFF, and ReHo images were co‐registered to the T1‐weighted structural image using boundary‐based registration and then normalized to a custom adolescent template using the top‐performing SyN registration provided by ANTs (Avants et al., 2011; Ciric et al., 2021; Greve & Fischl, 2009). Finally, a gray matter mask was generated as the intersection between a gray matter mask from T1‐weighted images with 90% coverage over all subjects and overall coverage masks from registered pcASL and rfMRI scans.

Methodology for estimating pIMCo

For each subject, we calculated voxel‐wise coupling between CBF, ALFF, and ReHo images to produce one pIMCo image per subject. The full pIMCo estimation pipeline is summarized in Figure 1. On average, pIMCo estimation for one subject took approximately 5 min.

FIGURE 1

Step‐by‐step diagram of pIMCo estimation pipeline as described in the methods. Coupling images are generated from intermodal images for each subject individually. pIMCo estimation is performed at each voxel location across subject‐specific images

First, we applied the gray matter mask to each of the three modalities. Then, within each masked modality, we globally scaled intensities to a mean of 0 and a variance of 1 across all voxels in the gray matter mask. This scaling is necessary because eigendecomposition is later performed on local covariance matrices; if modalities are defined on drastically different scales, decomposition outputs would reflect differences in scale between modalities rather than local covariance structures. Next, for each voxel, we extracted local neighborhoods from each of the three modalities and weighted voxels within these local neighborhoods proportional to a Gaussian kernel over their Euclidean distances from the central voxel – in our study, we used FWHM = 3 mm, which corresponds to 7 × 7 × 7 voxel (14 × 14 × 14 mm) local neighborhoods and a standard deviation of 1.27 mm for the Gaussian kernel, as prior WLR‐based IMCo studies have found this FWHM to be informative (Valcarcel, Linn, Vandekar, et al., 2018b). We also show that the FWHM can be varied parametrically in the Supplementary Materials; however, with even a small increase in FWHM to 5 mm, which corresponds to 11 × 11 × 11 voxel neighborhoods, pIMCo estimation for one subject took approximately 12 min.

Then, we calculated the 3 × 3 weighted covariance matrix between the neighborhoods, performed eigendecomposition on it using the eigen() function in R, and extracted the proportion of variance explained by the first/largest eigenvalue.

Finally, once all the voxel‐wise proportional first eigenvalues were extracted, we linearly shifted and scaled these values such that their theoretical range was [0, 1] and performed a logit transformation. While this monotonic transformation makes the coupling value more challenging to interpret, it appropriately emphasizes extreme values of coupling and changes the domain of coupling values from the bounded domain to an unbounded domain , where is the total number of modalities. The transformation results in a coupling value that more closely follows a normal distribution and is thus expected to have improved behavior with post‐hoc statistical analyses (Supplementary Figure S1).

This resulted in our voxel‐level pIMCo image for that subject. For any voxel in that image, a large value suggests that the voxel's local covariance matrix across modalities could be well‐summarized in a single dimension while a small value suggests multiple dimensions would be necessary to characterize the covariance structure.

For reference, in the three‐modality setting, coupling values of −2, 0, and 2 correspond to the first eigenvector explaining 41%, 67%, and 92% of the total variance in that local neighborhood respectively. In the two‐modality setting, coupling values of −2, 0, and 2 correspond to the first eigenvector explaining 56%, 75%, and 94% of the total variance in that local neighborhood respectively.

Voxel‐wise statistical analysis

We created descriptive coupling maps by taking the means and variances across all 803 subjects' pIMCo values at each voxel location in volumetric space. We then projected these mean and variance images to the cortical surface using PySurfer for visualization of spatial heterogeneity and cortical patterns (Waskom et al., 2020).

To investigate the biological relevance of pIMCo, we used linear regression at each voxel to explore whether coupling was associated with age or sex. In all linear regressions, we controlled for in‐scanner motion for both ASL and rfMRI scans. To account for multiple comparisons in these voxel‐level tests, we controlled the false discovery rate at 5% (Benjamini & Hochberg, 1995). Then, we created binary thresholded masks indicating which voxels displayed a significant effect for each of age and sex. For this and following analyses, we performed identical modeling of each of the three modalities individually to explore whether age and sex effects were present and corresponded to the observed associations with pIMCo.

Spin testing

To visualize the extent of voxels where coupling was associated with age and sex, we counted the proportion of voxels with statistically significant age or sex effects in each of the Yeo 7 functional networks on the cortex as well as in subcortical regions in the Automated Anatomical Labeling (AAL) atlas (Thomas Yeo et al., 2011; Tzourio‐Mazoyer et al., 2002).

Next, we tested whether the proportion of significant voxels in each functional network was enriched when compared to the proportion of significant voxels overall. Because there is an underlying spatial distribution of significant voxels, we used the spin test (Alexander‐Bloch et al., 2018). Briefly, the spin test is a permutation‐inspired testing procedure that rotates the FreeSurfer sphere randomly to create an underlying null distribution that preserves spatial patterns. The null hypothesis is that there is no spatial enrichment of the proportion of significant voxels in the specified functional network compared to across the cortex overall. In our study, we estimated the null distribution over 2000 permutations – for each permutation, we recorded the Jaccard similarity index between the thresholded p‐value map and each of the Yeo 7 networks. Finally, for each network, we calculated the p‐value as the proportion of null Jaccard similarity indices equal to or greater than the observed Jaccard similarity index.

Code and data availability

An R package for calculating pIMCo images is available at: https://github.com/hufengling/pIMCo. All code for analysis is available at: https://github.com/hufengling/IMCo_analyses. Software and packages used for pre‐processing, pIMCo calculation, or analysis are cited in the references – R Core Team (2021); Wickham (2019); Xie (2021a); Xie (2021b); Allaire et al. (2021); Xie (2016); Wickham et al. (2019); Muschelli (2021a); Avants (2020); Kandel et al. (2020); Müller (2020); Fischl (2012); Schäfer (2021); Jenkinson et al. (2012); Muschelli et al. (2015); Ren (2021); Muschelli (2021b); Yushkevich et al. (2006).

PNC data are publicly available in raw format at https://www.ncbi.nlm.nih.gov/projects/gap/cgi-bin/study.cgi?study_id=phs000607.v3.p2.

RESULTS

Coupling varies across the cortical surface and is increased in discrete regions, bilaterally

We calculated voxel‐wise mean and variance maps of coupling values to characterize spatial patterns in CBF‐ALFF‐ReHo coupling and visualized these on the cortical surface. Throughout the cortical surface, all voxels, on average, showed strong coupling, and voxels with stronger average coupling also tended to have higher variance between subjects (Pearson correlation = 0.69). The average voxel‐wise mean coupling value was 0.99 (sd = 0.37; range = [0.27, 3.30]). The average voxel‐wise variance was 0.91 (sd = 0.20; range = [0.45, 2.64]).

Visual comparison of voxel‐wise descriptive maps with the Desikan‐Killiany cortical atlas (Desikan et al., 2006) suggested that coupling is especially strong in the following regions, bilaterally: superior frontal gyrus, paracentral gyrus, caudal anterior cingulate, posterior cingulate, isthmus cingulate, pericalcarine, lateral occipital, and insula (Figure 2). Comparing to the Yeo 7 functional networks, areas of strong coupling are observed primarily in the frontoparietal (p = 0.0125) and default networks (p = 0.039; Figure 3).

FIGURE 2

Coupling values are spatially heterogeneous across the cortical surface. (a) Voxel‐wise means across subjects of cortical coupling values between cerebral blood flow (CBF), amplitude of low‐frequency fluctuations (ALFF), and regional homogeneity (ReHo). Larger values indicate stronger coupling. (b) Voxel‐wise variances across subjects of cortical coupling values between CBF, ALFF, and ReHo

FIGURE 3

Proportion of voxels in automated anatomical labeling subcortical regions and Yeo 7 cortical networks that showed significant coupling‐age and coupling‐sex associations when in‐scanner motion was included as a covariate (FDR corrected p < 0.05). Spin test was performed for all Yeo 7 networks; significant p‐values are reported (p < 0.05)

These voxel‐wise descriptive maps of coupling showed unique information when compared to the descriptive maps of each of the individual modalities (Supplementary Figure S2).

CBF‐ALFF‐ReHo coupling evolves with age throughout gray matter structures

Linear associations between strength of coupling and age were present in subcortical structures and cortical networks (Figure 3; corrected p < 0.05). Figure 4 shows an example of such an association between coupling and age as well as individual modalities and age at one voxel in the default network. In subcortical structures, age‐related changes in CBF‐ALFF‐ReHo coupling occurred primarily in the caudate and pallidum, though such changes were also common in the hippocampus, putamen, and thalamus.

FIGURE 4

Example of associations between individual modalities and age as well as associations between coupling and age at a single voxel in the default network. Each point represents the value at that voxel for one subject. Best‐fit lines from univariate linear regression are shown

In cortical networks, coupling and age associations were rare in all networks except the frontoparietal and default networks (Figure 5). These two networks were also the networks in which the average strength of coupling across subjects appeared to be highest (Figure 2). Spin testing between functional networks and age‐related changes in coupling showed enrichment of coupling and age associations in the frontoparietal (p = 0.013) and default networks (p = 0.039).

FIGURE 5

(a)Thresholded maps of voxels with significant coupling associations with age after FDR correction at 0.05. (b) Thresholded maps of voxels with significant coupling associations with sex after FDR correction at 0.05

Similar patterns were observed when FWHM was varied to either 2 or 5 mm instead of 3 mm (Supplementary Figures S3 and S4). Of interest, the overall proportions of voxels showing significant coupling‐age associations was highest at FWHM of 5 mm, intermediate at FWHM of 3 mm, and lowest at FWHM of 2 mm. This trend suggests that using a larger FWHM may result in coupling values that are more sensitive to group differences as pIMCo estimation is smoothed over a larger local neighborhood, though a trade‐off with spatial specificity and computational resources must be considered.

CBF‐ALFF‐ReHo coupling varies between males and females, primarily in subcortical regions

Associations between CBF‐ALFF‐ReHo coupling and sex were present primarily in the hippocampus and thalamus (Figure 3; corrected p < 0.05). Sex differences in coupling were rare in other subcortical structures and all functional networks – only 1% to 3% of voxels in these regions showed coupling and sex associations (Figure 5). Spin testing between functional networks and sex‐related changes in coupling revealed enrichment of coupling and sex associations in the frontoparietal network (p = 0.012), despite the small proportion of the frontoparietal network that exhibited coupling associations with sex. Again, similar results were observed when the FWHM parameter was varied (Supplementary Figures S3 and S4).

pIMCo provides a consistent estimator of local coupling compared to WLR‐based IMCo

In 2016, Vandekar et al. introduced a method to study IMCo relationships at the single voxel level based on local weighted linear regression (WLR) slopes (Vandekar et al., 2016). Because this method relies on estimating WLR slopes between modalities, it is inherently limited to the two‐modality setting, cannot account for statistical relationship between modalities, and requires specification of one modality as the independent variable, leading to asymmetry. These limitations are demonstrated by a two‐modality example in Figure 6. This figure was generated by choosing a random voxel with high coupling, as estimated by pIMCo, from a random subject and then calculating the WLR‐based IMCo value with ALFF as the independent variable (WLR‐ALFF), the WLR‐based IMCo value with CBF as the independent variable (WLR‐CBF), and the pIMCo coupling value.

FIGURE 6

Two‐modality example showing pIMCo results in comparison to weighted linear regression (WLR)‐based IMCo results. Blue line represents coupling value from WLR‐based IMCo if amplitude of low‐frequency fluctuations is defined as the independent variable (slope = 0.25). Orange line represents coupling value from WLR‐based IMCo if cerebral blood flow is defined as the independent variable (slope = 1.25). Black line and ellipse represent principal component analysis (PCA) results; no reference specification is needed (coupling value = 2.13). Larger dot sizes correspond to increased weights in the WLR and weighted PCA

We see that the WLR‐ALFF coupling value is 0.25, indicating little coupling, but the WLR‐CBF coupling value is 1.25, indicating five times as much coupling. Thus, when using WLR‐based IMCo, two possible coupling values exist at every voxel, and there are no guarantees that analyses will show comparable findings between the two. Additionally, these WLR coupling values only describe the trend of the relationship between ALFF and CBF, but do not account for the statistical strength of that relationship.

In contrast, the pIMCo coupling value is 2.13 and does not require specification of which modality is treated as the independent variable, leading to a symmetric and consistent definition of coupling. This pIMCo value does not describe the effect size of the relationship between ALFF and CBF and instead describes the strength and shape of the relationship – a high value of coupling suggests that the shape of data from that neighborhood looks like a long ellipsoid, while a low value of coupling suggests that the shape is more spherical.

We then repeated this comparison for 100 randomly sampled voxels across 50 randomly sampled subjects for a total of 5000 voxels. Again, for each voxel, we computed the WLR‐ALFF value, the WLR‐CBF value, and the pIMCo coupling value. Figure 7 shows the results of this random sampling comparison.

FIGURE 7

Comparison between pIMCo coupling values and weighted linear regression (WLR)‐based IMCo coupling values across 5000 randomly sampled voxels. (a) Relationship between WLR‐amplitude of low‐frequency fluctuations (ALFF) and WLR‐cerebral blood flow (CBF) coupling values showing a strong inverse relationship with instability when WLR‐ALFF approaches 0 and when WLR‐CBF approaches 0. (b) Relationship between pIMCo and WLR‐ALFF coupling values. Color of points corresponds to the magnitude of the WLR‐CBF coupling value, with lighter color indicating higher magnitude. (c) Relationship between pIMCo and maximum WLR magnitude showing strong correlation. (d) Pairwise non‐parametric Kendall correlations between pIMCo values, WLR‐ALFF values, and maximum WLR magnitudes showing the pIMCo and maximum WLR magnitude have the strongest pairwise correlation

Panel A shows a strong inverse relationship between the WLR‐ALFF and WLR‐CBF values such that, if the WLR‐ALFF value is very close to 0, suggesting no coupling, the WLR‐CBF value is unstable and is either very small or very large; the vice versa is also true. Panel B shows that when the pIMCo value is small, the WLR‐ALFF value is also small; however, when the pIMCo value is large, the WLR‐ALFF value can be small, large, or close to zero, suggesting either strong coupling or no coupling. Thus, Panel B shows that pIMCo and WLR‐ALFF estimates can lead to seemingly contradictory results. Coloring the points by their respective WLR‐CBF magnitudes reveals that this contradiction is due to the inverse relationship shown in Panel A – voxels with large pIMCo values and zero WLR‐ALFF values have large WLR‐CBF magnitudes. Similar trends can be observed when plotting pIMCo values against WLR‐ALFF magnitudes and WLR‐CBF magnitudes (Supplementary Figure S5).

In Panel C, we plot pIMCo values against the maximum magnitude of the WLR‐ALFF and WLR‐CBF values at each voxel and see strong correlation between the two. Voxels with large maximum WLR magnitudes can be interpreted as voxels that either WLR‐ALFF or WLR‐CBF identifies as having strong coupling. This is further supported by Panel D, which displays pairwise Kendall correlations for pIMCo values, WLR‐ALFF values, and maximum WLR magnitudes and shows pIMCo values and maximum WLR magnitudes are highly correlated.

Overall, these results illustrate the asymmetry of WLR‐based IMCo estimates in the two‐modality setting and show that apparent contradictions between pIMCo and WLR‐based IMCo estimates are resolved when the asymmetry of WLR‐based IMCo is taken into account. Thus, pIMCo provides a more consistent estimate of local coupling than WLR‐based IMCo even in the two‐modality setting.

DISCUSSION

As growing emphasis is placed on the acquisition of multi‐modal data, new methodologies are necessary to enable these analyses. In this manuscript, we introduce pIMCo, a generalized approach to estimating local IMCo that can be applied to two or more modalities, can be interpreted as a direct summary of local covariance matrices, and is symmetric (Figure 6). We show that pIMCo not only allows for coupling analyses in more than two modalities, but also is superior to WLR‐based IMCo in the two‐modality setting. Our method can be used with any combination of volumetric images to produce single‐subject, voxel‐resolution coupling images which can then be analyzed using standard techniques. We applied our proposed method to show significant coupling between cerebral blood flow, resting state fluctuations, and local connectivity throughout the brain. We then used voxel‐level analyses to characterize how coupling varies in neurodevelopment.

Coupling of cerebral blood flow, resting state fluctuations, and local connectivity vary with age and sex

Neurovascular coupling measured by CBF and ALFF has been previously thought of as a measure of the process by which neuronal activity induces changes in local vasculature to support nutrient demands (Armstead, 2016; Presa et al., 2020). This neurovascular coupling has been shown to evolve over neurodevelopment as well as in disease (Baller et al., 2022; Girouard & Iadecola, 2006; Hu et al., 2019; Jin et al., 2020). Here, we extend this concept of neurovascular coupling and use CBF‐ALFF‐ReHo coupling to explore not only the relationship between neuronal activity and blood flow, but also how this relationship is attenuated or enhanced by the additional consideration of the degree of local neuronal activity coordination. Thus, regions with strong CBF‐ALFF‐ReHo coupling may potentially reflect improved local optimization, where neuronal activity at a voxel induces, or is induced by, additional neuronal activity within its local network, and together influence local changes in blood flow.

We found that CBF‐ALFF‐ReHo coupling was spatially heterogeneous and varied with age and sex in neurodevelopment in both subcortical structures and functional networks. These findings, which uncover otherwise undetectable intermodal interactions, are unique to those from individual modality analyses.

We noticed that regions with higher CBF‐ALFF‐ReHo coupling across subjects also tended to have higher variance in coupling. This suggests that these regions may be biologically interesting in the context of CBF‐ALFF‐ReHo coupling and important regions for future exploration, since they appear to demonstrate pronounced differences in coupling phenotypes between subjects and could be associated with other variables of interest, such as clinical phenotypes.

Spin testing showed that the high proportion of coupling‐age associations in the frontoparietal and default networks were enriched when compared to the cortex overall. This suggests that in neurodevelopment, there is change in not only blood flow, resting state fluctuations, and local connectivity individually, but also in the strength of interaction among these features. These findings are consistent with and fortify the literature demonstrating the importance of frontoparietal and default networks as regions for change in neurodevelopment (Baller et al., 2022; Chai et al., 2017; Fair et al., 2008; Lin et al., 2019). Outside of neurodevelopment, our findings are consistent with previous work showing that coupling between CBF and functional connectivity strength is stronger in frontoparietal and default networks than regions outside these networks (Tak et al., 2015). In subcortical structures, high proportions of coupling associations with age seen in the caudate, pallidum, hippocampus, and thalamus suggest that modulation of vascular, resting state fluctuations, and local connectivity coupling may be necessary in the development of movement, memory, and fundamental brain activities.

In adolescent neurodevelopment, in response to hormonal changes and environmental exposures, myelinogenesis remains active, regional neurocircuitry is strengthened, and brain plasticity is supported by neuronal proliferation, rewiring, and dendritic pruning (Arain et al., 2013; Baum et al., 2020; Giedd et al., 1999; Yurgelun‐Todd, 2007). Additionally, cerebral blood flow has been shown to decrease rapidly in early adolescence (Satterthwaite, Shinohara, et al., 2014b). Together, these age‐related microstructural changes facilitate more efficient inter‐neuron communication in terms of both nutrient demand and local network efficiency, thus changing the overall coupling between blood flow, neuronal activity, and local connectivity (Bercury & Macklin, 2015).

High proportions of coupling associations with sex in the hippocampus and thalamus suggest that male–female differences in memory and related cognitive functions between males and females could be explained in part by the strength of regional brain metabolism as measured by cerebral blood flow, resting state fluctuations, and local connectivity coupling. In the cortex, the rarity of coupling‐sex associations suggests that this coupling may not play a role in explaining cortical sex‐based differences or that this relationship is more complex than our analyses could uncover. These cortical findings are of interest when compared to previous work showing a high proportion of the cortex had significant associations between sex and CBF‐ALFF coupling (Baller et al., 2022). Together, these studies demonstrate that three‐modality coupling identifies unique relationships when compared to two‐modality coupling.

Notably, despite low coupling‐sex signal in the frontoparietal network, spin testing showed enrichment of coupling associations with sex in this network. Since spin testing is a spatial permutation test that uses the coupling‐sex association thresholded p‐value map to generate data under the null, this significant finding is likely due to a combination of overall rare coupling‐sex associations in the cortex and the spatial distribution of these associations within the frontoparietal network. This finding highlights potential shortcomings of using the spin test for enrichment analysis. Additionally, statistical analysis of enrichment in subcortical structures is not yet possible, so more methods development is needed in this area.

Limitations and future directions

pIMCo is designed to summarize a complex local covariance structure, which necessarily leads to outputs that cannot fully characterize the intricacies of coupling. This is especially true in settings with more than two modalities – the covariance structure becomes even more complex, and more information is lost when it is summarized. As such, high coupling values can result from many different covariance features, and it is challenging to understand the basis of these high values.

The symmetric nature of pIMCo is conducive to more consistent interpretation when compared to slope‐based IMCo, since pIMCo does not depend on reference modality specification. However, while slope is a biologically intuitive measure that can be interpreted as capturing the directionality and effect size of a relationship, the generalizable nature of the pIMCo value makes it challenging to interpret biologically, especially when used on more than two modalities. Instead, it is most accurately interpreted in statistical terms – as a measure of the proportion of variance explained by the first eigenvector, as a measure of how well the local covariance structure could be summarized in one dimension, or as a measure of how ellipsoidal the neighborhood is instead of spherical. This is shown in Figure 6, where the WLR coupling value estimates that, on average, a 1 unit increase in CBF corresponds to a 1.25 unit increase in ALFF, while the pIMCo coupling value describes the ratio of the ellipse's major axis to minor axis and suggests that the statistical strength of the relationship between ALFF and CBF is strong.

Next, since PCA functions as a linear dimension reduction technique, it is most effective at summarizing data whose shape is roughly ellipsoid. However, it is unlikely that such an assumption holds for all voxels and all combinations of different modalities. For example, in a two‐dimensional neighborhood, it could be possible that there is a strong quadratic relationship between modalities, but the data is not well‐summarized by one eigenvalue and its corresponding eigenvector. In such cases, an IMCo technique based on manifold learning concepts could be a useful improvement that picks up on otherwise undetected intermodal relationships.

Finally, pIMCo is designed to estimate coupling in cross‐sectional multi‐modality data sets. However, there is also rich covariance information in longitudinal data sets with one or multiple modalities. Future implementations of pIMCo suitable for longitudinal data could investigate whether IMCo changes over time in neurodevelopment or disease, and if so, whether such changes may be useful in identifying biomarkers for clinical phenotypes.

CONCLUSION

pIMCo offers a novel perspective for summarizing the overall covariance structure between more than two modalities as well as a generalized, symmetric approach for describing coupling in the two‐modality setting. Here, we applied this method to the analysis of coupling between cerebral blood flow, resting state fluctuations, and local connectivity images. This analysis revealed patterns in neurodevelopment with respect to age and sex that differed from those present in any individual modality. As multi‐modal data becomes more common, we hope that pIMCo will serve as a tool for capturing complex intermodal relationships and enable more sophisticated descriptive analyses, improved prediction efforts, and novel methodological advances, among others.

AUTHOR CONTRIBUTION

Fengling Hu: Conceptualization, methodology, software, validation, formal analysis, investigation, writing – original draft, writing – review and editing, visualization. Sarah M. Weinstein: Methodology, software, formal analysis, investigation, writing – review and editing. Erica B. Baller: Investigation, resources, writing – review and editing. Alessandra M. Valcarcel: Conceptualization, methodology, software, investigation. Azeez Adebimpe: Software, investigation, resources, data curation. Armin Raznahan: Conceptualization, writing – review and editing. David R. Roalf: Conceptualization, writing – review and editing. Timothy E. Robert‐Fitzgerald: Methodology, software, data curation, visualization. Virgilio Gonzenbach: Methodology, software, visualization. Ruben C. Gur: Conceptualization, investigation, resources, funding acquisition. Raquel E. Gur: Conceptualization, investigation, resources, funding acquisition. Simon Vandekar: Conceptualization, software, formal analysis, investigation, writing – review and editing. John A. Detre: Conceptualization, writing – review and editing. Kristin A. Linn: Conceptualization, software, formal analysis, investigation, writing – review and editing, funding acquisition. Aaron Alexander‐Bloch: Methodology, software, formal analysis, investigation, resources. Theodore D. Satterthwaite: Conceptualization, validation, investigation, resources, data curation, writing – review and editing, supervision, project administration, funding acquisition. Russell T. Shinohara: Conceptualization, methodology, validation, investigation, resources, writing – review and editing, supervision, project administration, funding acquisition.

FUNDING INFORMATION

This study was supported by grants from the National Institute of Mental Health (R01MH112847 to TDS and RTS; R01MH120482 and R01MH113550 to TDS; R01MH123550 to RTS; R01MH107235 to RCG; 2T32MH019112‐29A1 to EBB; R01MH120174, R01MH119185, and R56AG066656 to DRR). FH was supported by NIH Medical Scientist Training Program T32 GM07170. SMW was supported by the National Science Foundation Graduate Research Fellowship Program (NSF GRFP). The PNC data used were funded by NIMH RC2 grants (MH089983 and MH089924 to REG). Funding sources were not involved in study design, data analysis, manuscript preparation, or submission decisions.

CONFLICT OF INTEREST

RTS receives consulting income from Octave Bioscience and compensation for reviewership duties from the American Medical Association. AMV is currently employed by Genentech. The authors report no conflicts of interest.

Figure S1 Example of change in distribution due to the proposed transformation from proportional first eigenvalue to coupling value. Colored lines represent density plots of proportional first eigenvalues or coupling values from five randomly selected subjects. Black lines represent the corresponding normal distributions – mean and standard deviation was calculated for both proportional first eigenvalues and coupling values across the five subjects.

Click here for additional data file.

Figure S2 A) Voxel‐wise means across subjects of cortical CBF. B) Voxel‐wise variances across subjects of cortical CBF. C) Voxel‐wise means across subjects of cortical ALFF. D) Voxel‐wise variances across subjects of cortical ALFF. E) Voxel‐wise means across subjects of cortical ReHo. F) Voxel‐wise variances across subjects of cortical ReHo.

Click here for additional data file.

Figure S3 Results from re‐analysis using FWHM of 2mm. Proportion of voxels in AAL subcortical regions and Yeo 7 cortical networks that showed significant coupling‐age and coupling‐sex associations when in‐scanner motion was included as a covariate (FDR corrected p < 0.05). Spin test was performed for all Yeo 7 networks; the two lowest p‐values are reported for each of age and sex. No p‐values were significant at alpha level of 0.05, but the lowest p‐values correspond to similar networks as analyses done with FWHM of 3mm and 5mm.

Click here for additional data file.

Figure S4 Results from re‐analysis using FWHM of 5mm. Proportion of voxels in AAL subcortical regions and Yeo 7 cortical networks that showed significant coupling‐age and coupling‐sex associations when in‐scanner motion was included as a covariate (FDR corrected p < 0.05). Spin test was performed for all Yeo 7 networks; the two lowest p‐values are reported for each of age and sex. These p‐values correspond to similar networks as analyses done with FWHM of 2mm and 3mm.

Click here for additional data file.

Figure S5 Comparison between pIMCo coupling values and WLR‐based IMCo coupling values across 5,000 randomly sampled voxels. A) Relationship between WLR‐ALFF and WLR‐CBF coupling value magnitudes showing a strong inverse relationship. B) Relationship between pIMCo and WLR‐ALFF magnitudes. Color of points corresponds to the magnitude of the WLR‐CBF coupling value, with lighter color indicating higher magnitude. C) Relationship between pIMCo and WLR‐CBF magnitudes. Color of points corresponds to the magnitude of the WLR‐ALFF coupling value, with lighter color indicating higher magnitude. D) Pairwise non‐parametric Kendall correlations between pIMCo values, WLR‐ALFF magnitudes, and WLR‐CBF magnitudes.

Click here for additional data file.