Covariance and Correlation Analysis of Resting State Functional Magnetic Resonance Imaging Data Acquired in a Clinical Trial of Mindfulness-Based Stress Reduction and Exercise in Older Individuals

We describe and apply novel methodology for whole-brain analysis of resting state fMRI functional connectivity data, combining conventional multi-channel Pearson correlation with covariance analysis. Unlike correlation, covariance analysis preserves signal amplitude information, which feature of fMRI time series may carry physiological significance. Additionally, we demonstrate that dimensionality reduction of the fMRI data offers several computational advantages including projection onto a space of manageable dimension, enabling linear operations on functional connectivity measures and exclusion of variance unrelated to resting state network structure. We show that group-averaged, dimensionality reduced, covariance and correlation matrices are related, to reasonable approximation, by a single scalar factor. We apply this methodology to the analysis of a large, resting state fMRI data set acquired in a prospective, controlled study of mindfulness training and exercise in older, sedentary participants at risk for developing cognitive decline. Results show marginally significant effects of both mindfulness training and exercise in both covariance and correlation measures of functional connectivity.


INTRODUCTION
Evaluation of resting state fMRI functional connectivity (FC) currently is dominated by two methods: Seed-based correlation (SBC) (Shehzad et al., 2009;Zuo et al., 2010) and spatial independent component analysis (sICA) (Beckmann et al., 2005). SBC conventionally is evaluated by Pearson correlation of two time series extracted either from seed regions of interest (ROIs) or voxels, essentially as first described by Biswal et al. (1997). Pearson correlation is a dimensionless, normalized measure that is invariant with respect to signal amplitude. sICA is intrinsically insensitive to signal amplitude. Accordingly, commonly used analysis procedures normalize input time series to unit variance as a preliminary step (e.g., Beckmann et al., 2005;Allen et al., 2011). Concurrently, substantial evidence has accumulated during the past decade indicating that the amplitude of spontaneous blood oxygen level dependent (BOLD) fluctuations is a meaningful indicator of the brain's functional integrity in psychiatric conditions (Gong et al., 2020), age-related cognitive decline (Vieira et al., 2020), and neurodegenerative diseases (Pan et al., 2017). Moreover, it has been reported that, in healthy individuals, the temporal standard deviation of the BOLD signal (SD BOLD ) indexes cognitive capacity in young as well as older individuals (Garrett et al., 2013;Grady and Garrett, 2014;Pan et al., 2017). This feature of the BOLD signal is ignored in conventional SBC and sICA.
Another issue addressed in this work concerns the dimensionality of the BOLD signal, which is known to be limited (Cordes and Nandy, 2006;Gotts et al., 2020). Functional connectivity (FC) is a second order statistic (Liegeois et al., 2017) that exists in a space of enormous dimension. Thus, given n ROIs, there are n · (n − 1)/2 unique ROI pairs and a corresponding number of potential FC measures. Thus, for example, if n = 300, as in the present work, n · (n − 1)/2 = 44,850. One strategy for dealing with this high dimensionality is to restrict the FC analysis to seed ROIs representing one or a very small number of a priori selected functional systems. This approach is suitable for testing a priori hypotheses concerning particular functional systems or loci within the brain. However, this option does not apply when the objective is to conduct a wholly data-driven, whole-brain FC study. Rational approaches to dealing with the discrepancy between the measured vs. true dimensionality of FC data have not been widely adopted.
Here, we present an approach to the challenge of obtaining whole-brain FC measures that incorporates dimensionality reduction while simultaneously accounting for the amplitude of spontaneous BOLD signal fluctuations. To this end we analyze resting state fMRI data acquired during the course of a large scale, prospective study of mindfulness meditation and exercise in older (age 65-84 years), sedentary participants at risk for developing cognitive decline. This study was conducted by the MEDEX (Mindfulness, EDucation, and EXercise) Research Group, a consortium comprising Washington University in Saint Louis (WUSM) and the University of California in San Diego (UCSD) and is the first study of its kind. It is registered in ClinicalTrials.gov (NCT02665481). A full description of the MEDEX study design is given in Wetherell et al. (2020). The rationale underlying the MEDEX study is that pharmacological treatments that halt or reverse aging-associated cognitive decline are not available. However, substantial evidence indicates that physical exercise ameliorates the manifold negative consequences of aging (Gallaway et al., 2017;Dauwan et al., 2021). Other studies have suggested that behavioral interventions, especially, mindfulness-based stress reduction (MBSR) (Kabat-Zinn et al., 1992) may improve cognitive function and reduce stress or depression in older individuals (Hazlett-Stevens et al., 2019). The question, then, is whether the effects of mindfulness training and exercise are detectable by analysis of resting state BOLD signals.
The methodology described herein was developed, in part, to address this question. We demonstrate novel methodology that addresses the problem of dimensionality in FC analysis while accounting for the amplitude of spontaneous BOLD signal fluctuations. We report resting state functional magnetic resonance imaging (rs-fMRI) outcomes derived from the MEDEX study. Non-neuroimaging outcomes of the MEDEX study will be reported elsewhere.

Participants and Study Design
Participants were recruited at two-sites, Washington University in Saint Louis (WUSM) and University of California San Diego (UCSD). Inclusion criteria included age 65-84 years, sedentary lifestyle, self-reported cognitive complaints but non-demented cognitive status, and no contraindication to magnetic resonance imaging (MRI), e.g., metal implants. Participants with biomarker evidence of preclinical Alzheimer's disease were not excluded. All participants gave written informed consent and received no remuneration. The IRB committees at WUSM and UCSD provided oversight over all aspects of the study.
Participants were randomly assigned to one of 4 interventions for an 18-month period, according to a 2 × 2 factorial design: (i) MBSR-only: Weekly instructor-led MBSR group-based classes for 10 weeks and then monthly booster sessions; (ii) Exercise-only: Twice-weekly instructor-led exercise group classes, including aerobic, strength, and functional training, for 6 months and then weekly booster session; (iii) MBSR + exercise: Both MBSR sessions and exercise sessions; (iv) Health education: instructorled sessions with health education content which included neither MBSR nor exercise and which matched the MBSR condition in session frequency and time. Participants were also instructed to practice at home over the entire 18-month duration of the study. The goal of home practice was up to 45 min daily mindfulness practice in the MBSR condition and 150 min/week exercise in that condition. MRI scanning was performed at baseline (before any intervention), at 6 months, and at 18 months.

Magnetic Resonance Imaging Acquisition
Two Siemens (Erlangen Germany) scanners equipped with 20channel head coils were used at WUSM. At UCSD, MRI was acquired using a GE MR750 3T MRI scanner (GE, Milwaukee, WI) equipped with an 8 Channel head coil (Table 1). Structural imaging included T1-weighted (WUSM MP-RAGE; TR = 2,400 ms, TE = 3.16 ms, TI = 1,000 ms; 1 × 1 × 1 mm voxels) (UCSD MPRAGE; TE = 3.036, TI = 1000 ms; and 1 × 1 × 1 mm voxels) and T2-weighted (WUSM SPACE; TR = 3,200 ms, TE = 458 ms; 1 × 1 × 1 mm voxels) (UCSD CUBE; TR = 2,500, TE = 73.37 ms; 1 × 1 × 1 mm voxels) anatomical images. Resting state fMRI (rs-fMRI) was acquired with a multi-echo sequence (WUSM TR = 2,960 ms, TE = 15, 31.3, 47.6, 63.9 ms; 4 × 4 × 4 mm voxels) (UCSD TR = 2,740 ms, TE = 14.8,28.4,42,55.6; 4 × 4 × 4 mm voxels) including 140 frames (volumes) per run. Up to 4 rs-fMRI runs were obtained in each session. During rs-fMRI acquisition, participants were shown a silent video of neutral content (relaxing nature scenes) and were instructed to keep their head still, stay awake, and not meditate. To simplify statistical comparisons, the present analysis includes only participants who completed all 4 resting state fMRI runs (23.3 min total WUSM, 25.6 min total UCSD) in all three scanning sessions. In accordance with the longitudinal experimental design, each participant was scanned with the same scanner during all three visits. All MRI scans were conducted at least 48 h after the participant's most recent exercise session (in class or at home) to avoid acute exercise effects in scan findings. Demographic information broken down by scanner is listed in Table 2. Cognitive performance and adherence data for the 4 intervention groups are listed in Table 1. Reported values are treatment group means. Cognitive Score is the normed (mean = 100, SD = 15), Fluid Cognition Composite test score from the NIH Toolbox Cognition Battery (Heaton et al., 2014), measured at baseline. Adherence is% classes attended over the study duration. MBSR and Health Education classes were once weekly for 10 weeks, then once monthly (total = 100 classes over 18 months). Exercise classes were twice weekly for 6 months, then once weekly (total = 25 over 18 months). Listed age refers to the baseline session. Ntot is participants scanned in all 3 sessions. Nret is participants contributing to the Results, retained after exclusion owing to excessive head motion in any session. %frames retained refers to data contributing to the Results after motion "scrubbing" (see below).

Resting State Functional Magnetic Resonance Imaging Processing
Pre-processing largely following methods described by Raut et al. (2019). Initial preprocessing was computed on data summed over all echoes and included rigid body correction for head motion within and across fMRI runs, correction of bias field inhomogeneities using the FAST module in FSL (Zhang et al., 2001), and slice timing correction. Atlas transformation was computed by 12-parameter affine registration of the structural T1w images to composition of affine transforms linking the fMRI data (head motion corrected functional frame average) to the atlas-representative target image (711-2B version of Talairach space). A scanner-specific target was generated for each of the three scanners (Buckner et al., 2004) to eliminate systematic atlas transform differences arising from variable T1w contrast. Transforms linking the functional data to the atlas representative target via the structural images were composed (fMRI→T2w→T1w→ atlas) and then applied in one-step to resample the functional data (4 echoes per frame) in atlas space (3 mm 3 voxels).
The atlas-transformed, multi-echo data were modeled according to standard theory (Poser et al., 2006) in which reconstructed image intensity depends mono-exponentially on TE. Thus, where S is intensity and S 0 is intensity extrapolated to TE = 0. S 0 and R * 2 are free parameters determined on the basis of multiple echoes (4 in this case). S 0 and R * 2 were estimated according to Eq. 1 separately for every voxel and frame using log-linear fitting. Empirical evidence (Power et al., 2018) shows that fluctuations in the value of S 0 primarily reflect spin history artifacts generated by head motion (Friston et al., 1996), whereas R * 2 reflects BOLD contrast (Ogawa et al., 1992) as well as changes in arterial pCO 2 (Birn et al., 2006). Accordingly, frame-to-frame variability in S 0 was eliminated by replacing time-dependent values with the (voxel-wise) fMRI run temporal average. The multi-echo modeling procedure then evaluated Eq. 1 at the TE corresponding to the second echo (31.3 ms for WUSM, 28.4 ms for UCSD) and output a volumetric time series that we here designate "Sfit." The Sfit volumetric time series acquired over 4 runs in each session were virtually concatenated. Next, to enable interpretation of fMRI signal fluctuation on an absolute scale, the functional data were intensity normalized (one scalar multiplier) to obtain a whole-brain mode value of 1,000. Thus, following mode-1,000 intensity normalization, a voxelwise temporal standard deviation of 10 corresponds to 1% rms signal modulation. Denoising began by marking frames for subsequent exclusion from the FC computations by reference to the DVARS timeseries, i.e., root-mean-square inter-frame intensity changes (Smyser et al., 2010;Power et al., 2014). The frame censoring criterion was adjusted on a per-session basis to accommodate baseline shifts in the DVARS measure (White et al., 2020). Frame censoring statistics are included in Table 2. The concatenated data then were spatially filtered (6 mm FWHM in each cardinal direction) and temporally filtered (demeaned, detrended, low-pass cut-off at 0.1 Hz). Additional denoising was accomplished by regression of timeseries using a strategy similar to CompCor (Behzadi et al., 2007). Imagederived nuisance regressors were extracted from FreeSurfer 6.0.0-segmented regions (Fischl, 2012) following co-registration with the functional data in atlas space. Nuisance regressors included (i) six rigid body parameter time series derived from within-run head motion correction; (ii) image-derived regressors extracted from multiple sub-regions within three anatomical compartments: white matter, ventricles, and the extracranial cerebrospinal fluid (CSF); (iii) the global signal averaged over the whole brain (Fox et al., 2009;Power et al., 2017;Ciric et al., 2018). Image-derived nuisance regressors were dimensionality reduced by PCA as previously described (Raut et al., 2019). The final number of nuisance regressors used to denoise the Sfit data varied according to the quality of the data (mean = 38, SD = 13, max = 85, min = 9). To account for the effects of the video stimuli shown during fMRI acquisition, the mean sessionspecific video fMRI response averaged over all participants was subtracted from each participant's timeseries. The preprocessed and denoised timeseries was extracted from 300 functionally defined brain regions of interest (ROIs) (Seitzman et al., 2020; Figure 1), and pairwise region of Interest (ROI) correlation values were computed, omitting frames previously marked for censoring.

Dimensionality Reduction and Derivation of Fixed Basis
The dimensionality of whole brain resting state fMRI data (number independent signals distinguishable from noise) is limited (Cordes and Nandy, 2006;Gotts et al., 2020). Accordingly, it is possible that dimensionality reduction may enhance sensitivity to experimental interventions by removing variance unrelated to resting state network (RSN) structure. Here, dimensionality reduction was effected by proper orthogonal decomposition (POD), an analytic technique closely related to principal component analysis (PCA) (Liang et al., 2002). We refer to this method as POD to emphasize that the objective is to derive a basis of reduced dimensionality on which to represent a high dimensional process. We retained the top 20 components, which is comparable to the number of non-noise components identified in prior work (Allen et al., 2011). Let X i represent the fMRI data where i indexes a particular session of a particular participant. X i is m × L i , where m is the number of ROIs (300 in the present work) and L i is number of rs-fMRI samples (total length of resting state data excluding censored frames) in session i. The covariance matrix of X i is The mean covariance matrix in the studied cohort isC = (1/N) where the eigenvectors of C constitute the columns of W and is a diagonal matrix of eigenvalues. The dimensionality reduced mean covariance matrix,Ĉ, is obtained by truncating , retaining the left upper 20 × 20 submatrix,ˆ . Thus, Thus,Ŵ constitutes a fixed basis of reduced dimensionality that provides a means of representing the covariance structure of all participants in a canonical format.
Left multiplying X i byŴ T yieldsŶ i =Ŵ T X i , the projection of session i's fMRI data onto the fixed basis. The covariance matrix of this projection is where C i is the full covariance matrix of session i. POD of this matrix yields In practice, this identity is only approximate (Ŵ i ≈Ŵ). Nevertheless, we define the diagonal entries of (1/L i )Ŷ iŶ T i ≡ˆ i as the estimated magnitude of 20 covariance components in session i and definê C i ≡Ŵˆ iŴ T as the projection of session i's covariance structure onto the fixed basis. The extent to which the diag (1/L i )Ŷ iŶ T i is equivalent to the first 20 eigenvalues of C i is the extent to whicĥ W TŴ i is equal to the 20 × 20 identity matrix, I. Similarity of subgroup eigenstructure is reported in Supplementary Figure 1.
Importantly,ˆ i andĈ i =Ŵˆ iŴ T are informationally equivalent (becauseŴ is fixed). Hence, ˆ i can be subjected to algebraic operations, e.g., averaging over participant subgroups and inputting into linear regressions. The results of these operations are linear combinations of ˆ i which can be inserted into the form of Eq. 4 and displayed as covariance matrices.
For notational simplicity, define i = diag ˆ i T . Thus, i represents the 20 covariance components corresponding to a particular visit of a particular participant reshaped as a 1 × 20 row vector. Let * denote some linear combination of { i }; then there exists a one-to-one correspondence between * and C * =Ŵˆ * Ŵ T . The asterisk in the preceding expression denotes any particular subgroup (e.g., participants in the exercise group imaged at 6 months). We present results using bothĈ * and * representations. To equalize total BOLD power over scanners, the individual covariance matrices were scaled by a site-specific factor ensuring that the traces of site-specific mean covariance matrices were equal (see Supplementary Figure 2). Further, sitespecific contributions to i were removed by linear regression.
Dimensionality reduction and basis definition for correlation (as opposed to covariance) FC are essentially similar. Thus, letr be the mean correlation matrix averaged over all participants and sessions. Thenr =ŵλŵ T is the dimensionality reduced mean correlation matrix and w contains eigenvectors defining the correlation FC basis. Equation 6 is analogous to Eq. 4. The remainder of the abovediscussed considerations, in particular, projection of individual FC components onto a fixed basis (Eq. 5), apply equally well to FIGURE 1 | ROIs projected onto the cortical surface. The total number of ROIs is 300. Each ROI is associated with one of 16 RSNs. "Unassigned" refers to regions of the brain in which fMRI signals are unreliable owing to susceptibility dropouts.
the case of correlation FC. Table 3 lists all variables used in the statistical evaluation and presentation of results.

Statistical Testing
Following dimensionality reduction, the covariance and correlation measures obtained in a particular session are represented as i and P i , respectively. To define additional nomenclature by example, let visit2−1 represent a longitudinal change operator. Thus, i.e., the effect of mindfulness training on longitudinal change in covariance FC over the first 6 months, and similar expressions We took the L 1 norm (sum of absolute values) of these 1 × 20 quantities as the measure of interest. Statistical significance of quantities of the form represented by Eq. 7 was assessed by permutation resampling over 10,000 repetitions. Thus, the link between participant and treatment was randomly shuffled, maintaining constant treatment group sizes, and the distribution of L 1 norm values compiled over repetitions. The likelihood of observing the true experimental outcome by chance then corresponds to a particular percentile of the surrogate distribution.

Covariance:Correlation Matrix Global Proportionality
As will shortly be shown,Ĉ andr evaluated over the full dataset exhibit strikingly similar "matrix topographies" (see section "Results"), i.e., differ, to a good approximation, by only a scalar factor. This relation can be symbolically represented asĈ ≈ υ ·r. where υ is a scalar factor. We fit this model by minimizing error over matrix block averages. Thus, we minimize where subscript k indexes matrix blocks and the bracket notation denotes averaging over entries within block. These blocks (delineated by heavy lines in Figure 2) are square on the diagonal and rectangular off the diagonal. The "blockwise average" approach to evaluating υ follows from the demonstration that dimensionally reduced FC matrices retain almost all RSN structure. The ordinary least squares estimate for υ is k< C> k <r> k / k <r > 2 k . The

Pearson correlation betweenĈ andr block averages is
Frontiers in Neuroscience | www.frontiersin.org Pearson correlation to avoid overloading the symbol r, which refers to fMRI signal correlations). As is true of Pearson correlation generally, the fraction of total variance accounted for by the global covariance:correlation proportionality model is η 2 . Theoretically, the value of υ depends sequence details as some sequences could weight white matter, gray matter and CSF differently. However, the impact of such dependencies likely are minor as the present methodology includes variance equalization across scanners (see Supplementary Figure 2). Figure 2 shows the whole-cohort, mean covariance and correlation matrices before (C,r) and after (Ĉ,r) projection onto their respective fixed bases. The block structure of the matrices, shown in Figure 2 replicates established findings reported in multiple rs-fMRI studies (Laumann et al., 2015;Gotts et al., 2020;Seitzman et al., 2020). It is visually evident that the block structure of these matrices is nearly unaffected by dimensionality reduction. The squared Pearson correlation betweenC k andĈ k is 0.985. The squared Pearson correlation betweenr k andr k is 0.996. Thus, dimensionality reduction preserves the block structure of both covariance and correlation FC matrices. At the same time, projection accounts for only 38.7% of total variance inC and 36.5% of total variance inr. Thus, variance outside the fixed bases resides almost entirely within subcomponents of RSNs. Figure 3 shows scree plots and eigenvectors (ROI weights) corresponding toĈ andr. The ROI weights in the basis vectors reflect the organization of major functional systems. Thus, the first component of bothŴ andŵ is dominated by the DMN and the second component is dominated by somatomotor dorsal, somatomotor ventral, and cingulo-opercular networks (SM_Dors, SM_Vent, Co.). Higher components exhibit progressively less RSN structure. This loss of RSN structure, coupled with asymptotically small eigenvalues as component indices approach 20, implies that the dimensionality reduction largely preserves meaningful variance in the original data. This point is addressed also in Supplementary Figure 1, which shows that the fixed bases reasonably well represent the correlation structure of participant subgroups.

Anatomical Topography of Covariance and Correlation Basis Vectors
Global Covariance: Correlation Proportionality Figure 4 shows the dimensionality reduced covariance and correlation matrices corresponding to the whole study cohort. It is evident that these matrices exhibit strikingly similar "matrix FIGURE 3 |Ĉ andr bases (eigenvectors) obtained by POD of the cohort-mean covariance and correlation FC matrices. Scree plots are shown above corresponding eigenvectors. The correlation eigenvectors have been multiplied by the scalar constant (υ = 3.67) that minimizes error in the global, block-wise proportionality model, C k = υ ·r k . The first few eigenvectors exhibit clustered, large magnitude loadings within resting state networks, e.g., the DMN. Network structure becomes fragmented at higher eigenvector indices in both the covariance and correlation representations of FC.
topographies." Moreover, the covariance-FC and correlation-FC basis vectors are similar (see Supplementary Figure 3). Blockwise fitting the model,Ĉ = υ ·r + , yielded a global covariance:correlation ratio (υ) of 3.67. The proportion of modelconsistent variance (η 2 ) is 0.866. Figure 4C shows the difference, C − υ ·r, i.e., focal deviations from global covariance:correlation proportionality. Such deviations include somewhat greater covariance:correlation ratios in parts of Visual network and the DMN and somewhat lesser ratios in somatomotor and cinguloopercular cortex. Although these deviations are potentially of physiological interest, they are quantitatively minor (13%).
Hence, further consideration of these focal deviations from global proportionality is deferred to future work.

Effects of MEDEX Interventions on Covariance and Correlation Functional Connectivity
The MEDEX study design includes two interventions (MBSR vs. no MBSR) × (Exercise vs. no Exercise) and three visits (Baseline, 6 months, 18 months), which generates 6 potential Intervention × Visit contrasts. We elected to simplify the analysis by focusing on first-order effects of the two types of intervention. Thus, study group comparisons were (1) exercise vs. no exercise (i.e., the exercise-only and MBSR + exercise groups, vs. the MBSR-only and Health Education groups), and (2) MBSR vs. no MBSR (i.e., the MBSR-only and MBSR + exercise groups, vs. the exercise-only and Health Education groups). Display of both covariance and correlation matrix results corresponding to all potential contrasts is not feasible. However, matrix results and component magnitude differences corresponding to contrasts yielding statistically significant results (p < 0.05, uncorrected) are shown in Figures 5, 6. The matrix and component magnitude displays are arranged in a 3 × 3 array with contrast over time in columns and contrast over intervention in rows. Figure 5 shows the effect of mindfulness training on covariance FC change at 18 months vs. baseline. Figure 6 shows the effect of exercise on correlation FC change at 18 months vs. baseline. Permutation testing of these effects is illustrated in Figure 7. Summary statistics covering all 6 condition contrasts for both covariance and correlation FC are listed in Table 4.

DISCUSSION
We report FC results obtained with novel methods combining conventional Pearson correlation FC with methodology that preserves BOLD amplitude information. This is the first large scale, prospective study of the effects of mindfulness training and exercise on resting state BOLD fMRI. In the present data, some effects of these interventions are formally significant (p < 0.05), omitting correction for multiple comparisons (Table 4). With correction for multiple comparisons (6 tests), these results are properly viewed as marginally significant. Several prior resting state fMRI studies concerned with the effects of MBSR or similar practices have reported increased conventional (Pearson correlation) FC between dorsolateral prefrontal cortex (dlPFC) and an array of other regions of the brain (Creswell et al., 2016;Taren et al., 2017;Kral et al., 2019). Similarly, exercise or aerobic fitness has been associated with increased (conventionally assessed) FC affecting a variety of widely distributed region pairs (Voss et al., 2010;McGregor et al., 2018). We did not specifically attempt to replicate those findings. However, the present evidence (Figure 6 and Table 4) suggests that the effects of exercise on conventional Pearson correction FC are exceedingly subtle and largely confined to areas of the cerebral cortex concerned with vision. Failure to replicate previously reported effects of exercise on FC measures has been explicitly noted before (Flodin et al., 2017).
Small participant samples undoubtedly account, in part, for replication failure. Sample sizes in all of the above-cited reports, except (Flodin et al., 2017), were a full order of magnitude smaller than the present one (Table 2). However, we suggest that principal challenge in data-directed functional connectivity studies is not fundamentally a matter of sample size. The fundamental problem is the dimensionality of the data space, which encompasses pairs of regions, numbering on the order of 10 4 in whole-brain studies utilizing dense spatial coverage (see Introduction). In practice, very little prevents an investigator focusing on selected region pairs, after which significant findings may emerge even after seemingly appropriate multiple comparisons correction. Dimensionality reduction offers a means of projecting whole-brain FC measures onto a space of manageable dimension. The investigator may vary the number components retained in the analysis but this does not bias the results provided that the eliminated components exhibit little evidence of structure (see Figure 3). Importantly, the dimensionality of resting state fMRI data is considerably smaller than the space of all densely sampled ROI pairs. This point has been made before (Cordes and Nandy, 2006;Gotts et al., 2020); it is here demonstrated in Figures 2, 3.
The present approach to the representation of FC by projection onto a fixed basis represents a greatly simplified version of previously published methodology (Madsen et al., 2017). In the present data, less than half (38.7%) of all BOLD variance is structured according to RSNs. This means that more than half is unstructured. Much of this unstructured variance arises from electronic noise (Liu, 2016). This variance enters into evaluations of seed-based FC but is not organized at the systems-level, hence, depresses measured correlations. The information content of the matrix and bar plot displays is identical. At 18 months vs. baseline, all participants, with or without MBSR training, showed increased activity in visual cortex and the DMN with more pronounced VIS-DMN negative covariance (upper matrices in right column). The 18-month change specific to MBSR was less focal and quantitatively modest (lower matrix in right column). In the notation of Eq. 7, the lower right bar plot represents MBSR−noMBSR vist3−1 . Statistical testing of | MBSR−noMBSR vist3−1 | is illustrated in Figure 7A (p = 0.03).
FIGURE 6 | 18-months vs. baseline effect of exercise on correlation FC. The format of this figure is identical to Figure 5. As in Figure 6, the information content of the matrix and bar plot displays is identical. Statistical testing of the lower right bar plot is illustrated in Figure 7B.
Dimensionality reduction eliminates variance not organized at the systems level, hence, avoids this problem. Additionally, projection of covariance and correlation measures onto a fixed basis provides a straightforward means of regressing out unwanted sources of variance, e.g., scanner dependencies. Parallel evaluation of covariance FC together with correlation FC is motivated by substantial evidence indicating that the amplitude of spontaneous BOLD fluctuations indexes cognitive capacity (Grady and Garrett, 2014;Pan et al., 2017). For example, it has been shown that BOLD signal variability is a correlate of age-related cognitive decline (Vieira et al., 2020). Moreover, it has been reported that the amplitude of spontaneous BOLD signal fluctuations correlates with performance measures, independently in young as well as older individuals (Garrett et al., 2013). Prior relevant work suggests that the amplitude of low frequency fluctuations (ALFF) is altered in adult long-term meditators (Berkovich-Ohana et al., 2016). More specifically, it has been reported that mindfulness or "mind-body" training  Table 4. | | and | P| denote the L 1 norm of quantities defined in Eq. 7. Squares roots were evaluated as this modification yielded a more normal (less right-skewed) distribution of permutation resampling surrogate values (see Figure 7). Red and blue highlight indicate p-values significant at p < 0.05 and p < 0.10, respectively, not corrected for multiple comparisons.
decreases ALFF in the default mode network (DMN) (Berkovich-Ohana et al., 2016;Yang et al., 2019) or hippocampus (Tao et al., 2019). [The hippocampus is closely linked to the DMN (Vincent et al., 2007)]. Such decreases in ALFF are broadly consistent with the results shown in Figure 5. Thus, decreased ALFF may not necessarily be inconsistent with a positive influence of MBSR on performance measures or mood. At the same time, focal ALFF increases, particularly in networks associated with cognitive control, have been reported as a correlate of mindfulness training (Wei et al., 2017), which is consistent with the work of Grady and Garrett (2014).
The quantity previously reported as SD BOLD (Garrett et al., 2013) or ALFF (Wei et al., 2017) is the square root of the quantity appearing on the diagonal of covariance FC matrices. Thus, covariance FC effectively comprises SD BOLD /ALFF but also broadens the analysis to include cross-RSN interactions which appear in off-diagonal blocks. Other laboratories have reported BOLD fMRI signal covariance matrices (Varoquaux et al., 2010). However, computation of BOLD fMRI covariance FC is hardly at all represented in the extant literature (but see Smyser et al., 2016). The present demonstration of approximately uniform global covariance:correlation proportionality (υ; Figure 4) is novel. This result implies that the amplitude of spontaneous BOLD fluctuations is approximately uniform over the cortical surface. Although the impact of MBSR on covariance FC was modest in the present data, Figure 5 shows covariance component magnitude changed over the 18 months between visits 1 and 3, unrelated to treatment. Specifically, BOLD fluctuations increased primarily in visual cortex and, to a lesser extent, the DMN, with concomitant strengthening of VIS:DMN negatively signed covariance. Much of this change is attributable to an increase in the magnitude of the first covariance component (see Figure 3). Conceivably, this change may reflect different videos played at visits 1 and 3, although the movie-evoked response was removed from the BOLD data before any FC computations. An explanation for this finding is not readily apparent. However, it does suggest that covariance FC (as well SD BOLD /ALFF) may be subject to variability owing to as yet poorly understood factors.

CONCLUSION
The principal objectives of this work are demonstration of the computational advantages of dimensionality reduction in the analysis of resting state fMRI data and introduction of parallel covariance/correlation analysis of functional connectivity. These advantages include (i) projection of FC measures, which, in "raw" form, exist in a space of very large dimension, onto a space of manageable dimension; (ii) enabling exclusion of unwanted sources of variance (e.g., scanner dependencies) using simple linear regression; (iii) exclusion from FC computations of unstructured variance that otherwise would bias FC measures toward zero; (iv) introduction of the global covariance:correlation proportionality constant, υ. Although the value of υ was unaffected by any of the present experimental manipulations, this ratio potentially carries physiological significance and may be revealing in other experimental contexts. Application of the novel methodology to the present data set revealed only marginally significant effects of both mindfulness training and exercise, in contrast to prior reports. In view of the unprecedentedly large participant sample relative to related prior work, this outcome raises questions concerning the replicability of prior findings.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Washington University School of Medicine IRB University of California, San Diego IRB. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
AS designed the study, developed analysis tools, performed the analysis, and wrote the manuscript. TN performed the analysis, created all figures, and edited the manuscript. JS designed the study, developed analysis tools, and edited the manuscript. EL designed the study, secured funding, and edited the manuscript. JW designed the study and secured funding. MV organized the study working group. JM provided statistical expertise. MY performed statistical analyses. DM provided data archiving. JG provided essential data transfer expertise. JR analyzed the data. DS oversaw data acquisition at the San Diego site. DB designed the study and secured funding. LE oversaw and performed data acquisition at the San Diego site. All authors contributed to the article and approved the submitted version.

FUNDING
This study was funded by NIH R01AG049689 by the National Institute of Aging together with the National Center for Complementary and Integrative Health, Office of Behavioral and Social Science Research, and the McKnight Brain Research Foundation. Additional support came from UL1TR002345 from the National Center for Advancing Translational Sciences. AS was supported by NIH U19 AG032438, R01 AG072694-01A1, and 1P30NS098577. DB was supported by R01 MH090786 from the National Institues of Health. JS and AS were also supported by P50 HD103525 to the Intellectual and Developmental Disabilities Research Center at Washington University.