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

Abstract 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.

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.

K E Y W O R D S
ASL, connectivity, coupling, fMRI, intermodal, MRI, neurodevelopment

| 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 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 higherlevel 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 lowfrequency 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 . 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 threemodality 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 locallycoordinated neuronal activity, which can be thought of as a more comprehensive measure of neurovascular coupling.

| 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 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 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: T R = 1810 ms; T E = 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: T R = 4000 ms; T E = 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: T R = 3000 ms; T E = 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 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 . Finally, BOLDweighted 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 timedomain 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.
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 voxelin 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. 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 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 wellsummarized 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 permutationinspired 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 permutationsfor 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) 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.
In cortical networks, coupling and age associations were rare in all networks except the frontoparietal and default networks F I G U R E 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) 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.

F I G U R E 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 F I G U R E 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

| 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; cor-

| 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.
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 relationshipa 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 S5). 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 singlesubject, 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.
4.1 | 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 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 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 . 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 modalitiesthe 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 termsas 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.

| 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 twomodality 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 inter-