Delayed cerebrovascular reactivity in individuals with spinal cord injury in the right inferior parietal lobe: a breath-hold functional near-infrared spectroscopy study

2.1. Participants

All participants were recruited from the New Jersey/New York metropolitan area. A total of 20 participants with SCI (14M, 6F, mean age = 46.3 ± 10.2 years) and 25 age– and sex-matched able-bodied (AB) controls (19M, 6F, mean age = 43.2 ± 12.28 years) participated in this study. Institutional Review Board (IRB) approval was obtained, and all participants provided written informed consent to the study. Participants were excluded from the study if any of the following applied to them: (1) within one-year post-spinal cord injury or in the acute phase of injury, (2) presence of tetraplegia or inability to perform upper limb motor movements, (3) history of or concurrent traumatic brain injury (TBI), (4) history of or concurrently has psychiatric disorders such as post-traumatic stress disorder, addiction, bipolar disorder, or schizophrenia, (5) presence of acute illness or infection, (6) history of chronic hypertension, diabetes mellitus, stroke, epilepsy or seizure disorders, multiple sclerosis, Parkinson’s disease, (7) illicit drug abuse within the past 6 months, (8) has Alzheimer’s disease or dementia, (9) not able to speak English, and (10) younger than 18 years of age or older than 65 years of age (D. Y. Chen, Di, Amaya, et al., 2024). The level of spinal injury, completeness of injury, and the duration since injury were collected from participants with SCI (Table 1).

2.2. Breath-hold Task Paradigm

All participants performed a breath-hold task paradigm consisting of 6 blocks of alternating periods of 15 s breath-hold and 30 s resting periods (Figure 1A). The task started off with 30 s of rest initially, during which the participants were instructed to breathe at their own comfortable pace in a relaxed manner. All instructions were displayed via a computer screen using E-Prime 2.0 software (Psychology Software Tools, Pittsburgh, PA), instructing the participants when to hold their breath and when to breathe normally. The words “rest” and “hold” were displayed in white font at the center of the screen on a black background. The participants also performed a resting-state task and N-back working memory task as part of a larger study; however, due to the focus of this study being cerebrovascular reactivity, only the breath-hold task results are analyzed and discussed here.

• Download figure
• Open in new tab

Figure 1.

A) fNIRS Breath-hold (BH) task paradigm consisting of 6 breath-hold sessions, 15s each, interleaved with 30s of rest. B) fNIRS Optode layout with 10 optode sources and 24 detectors configured across the medial prefrontal cortex, dorsolateral prefrontal cortex, motor cortex, and inferior parietal lobe. A total of 26 channels are used with a source-detector distance of 30mm and 4 short separation channels are used with a source-detector distance of 8.4mm, as denoted by S1, S2, S3, and S4.

2.3. fNIRS Data Acquisition

A continuous wave fNIRS system was used for all fNIRS data acquisition at a sampling frequency of 25 Hz, using both 690 nm and 830 nm of light (CW6 System, TechEn Inc., Milford, MA). A total of 10 optode sources and 24 detectors were arranged across the head to construct 26 long source-detector separation channels and 4 short channels (Figure 1B). The long channels were placed at a source to detector distance of 30 mm while the short channels were placed at a source to detector distance of 8.4 mm (Brigadoi & Cooper, 2015). These short channels were placed bilaterally on the left and right inferior parietal lobes and in the medial prefrontal cortex to reduce the impact of physiological noise from extracerebral regions on the data of interest. The 26 long channels covered the following regions of interest (ROI): the left and right dorsolateral prefrontal cortex (LDLPFC, RDLPFC, channels 1-4 and channels 5-8, respectively), medial prefrontal cortex (MPFC, channels 9&10), left and right motor cortex (LMOTOR, RMOTOR, channels 19-22 and channels 23-26, respectively), and left and right inferior parietal lobes (LIPL, RIPL, channels 11, 12, 15, & 16, and channels 13, 14, 17, & 18, respectively). The ROI were identified on each participant’s head using a neural navigation system (Brain Sight Neural Navigator, Rogue Research Inc., Canada) with the Montreal Neurological Institute’s (MNI-152) average structural brain template as a reference.

2.4. fNIRS Data Preprocessing

The fNIRS data were preprocessed in the following order: (1) raw light intensity values were converted to optical density values, (2) wavelet-based head-motion correction was performed using a threshold of 1.5 times the interquartile range (Molavi & Dumont, 2012), (3) band-pass filtering was applied using a pass-band from 0.01 to 0.15 Hz, (4) data were then converted to oxygenated hemoglobin (HbO), deoxygenated hemoglobin (HbR), and total hemoglobin (HbT) using the modified Beer-Lambert law (Delpy et al., 1988), (5) short channel regression was applied, and (6) channels were excluded based on the relative power in the task frequency range (Figure 2). Short channel regression was implemented using the proximal channel method, in which a short channel regressor was chosen for each long channel based on its proximity to the long channel (Gagnon et al., 2012). In this case, channel S1 was used to regress out noise from channels 1-4, & 9, channel S2 was used to regress out noise from channels 5-8, & 11, channel S3 was used to regress out noise from channels 11, 12, 15, 16, & 19-22, and channel S4 was used to regress out noise from channels 13, 14, 17, 18, & 23-26 (Figure 2). The short channel regression models were solved using the ordinary least squares method and were formatted as

• Download figure
• Open in new tab

Figure 2. fNIRS preprocessing pipeline applied to both SCI and AB groups.

The preprocessing steps were carried out in the following order: 1) convert raw intensity data to optical density data, 2) perform wavelet-based head motion correction with 1.5 times the interquartile range threshold, 3) band-pass filter from 0.01 to 0.15 Hz, 4) convert the filtered optical density data to oxy-hemoglobin (HbO), deoxy-hemoglobin (HbR), and total-hemoglobin (HbT) data, 5) regress the short source-detector channel from the long channels, and 6) perform channel exclusion using the relative power in the breath-hold frequency range. The short source-detector channel regression and channel exclusion steps were applied to all hemoglobin species.

where Y_{long channel} represents the fNIRS data from one of the long channels, X_{short channel} represents the fNIRS data from the proximal short channel, and ε represents the residuals. Short channel regression was performed on all subject data, across each channel separately and the residuals were kept for further analyses (Eq. 1).

After short channel regression, fNIRS channels were excluded based on the relative power in the breath-hold frequency range. Since the breath-hold task period occurred every 45 s (15 s breath-hold and 30 s rest), the dominant breath-hold task frequency range was defined as 40 s to 65 s, or 0.0154 Hz to 0.025 Hz. Therefore, the relative power was computed as the average power in the breath-hold frequency range (0.0154 to 0.025 Hz) divided by the total average power in the frequency range from 0 to 12 Hz (up to the Nyquist frequency) (Eq. 2). The average power was computed using the bandpower function in MATLAB 2024a (The MathWorks Inc., Natick, Massachusetts). Channels with less than 50% relative power in the breath-hold frequency range were considered “insufficient” and were excluded from subsequent analyses (Zvolanek et al., 2023). This method of channel exclusion not only allows us to remove noisy data, but also offers a measure of task compliance, since those who did not comply with the task would likely not show high power in the breath-hold frequency range. Short channel regression and channel exclusion were applied to HbO, HbR, and HbT separately, and sufficient data from each hemoglobin species were analyzed in the subsequent steps. Data from the start of the breath-hold task (beginning of the first 15s breath-hold period) to the end of the entire task (including the last 30s rest period) were analyzed. The first 30s of rest were discarded since it was included for the participants to relax and serve as a buffer period before the breath-hold stimuli began. All data preprocessing scripts were written and executed in MATLAB 2024a (The MathWorks Inc., Natick, Massachusetts) and used HOMER3 functions (Huppert et al., 2009).

2.5. Principal Component Analysis

After the fNIRS data were preprocessed, all data were standardized to z-scores by subtracting each timepoint by the mean of the entire time-series and then dividing by the standard deviation. Z-score standardization was performed for all channels and all subjects, for each hemoglobin species concentration separately. Principal component analysis (PCA) using a singular value decomposition algorithm was then applied to the z-scored data for each channel, in which the data from subjects in both AB and SCI groups were aggregated together in a [timepoint x subject] matrix. The square root of the eigenvalues of the covariance matrix was taken and multiplied by the eigenvector to obtain the principal component loadings. Since the first principal component explained the most variance, the coefficients of the first principal component (PC1 loading) were analyzed further. The PC1 loading value was inferred to represent the breath-hold response and serves as a measure of CVR amplitude, in which a higher PC1 loading represents a greater CVR amplitude response. PCA is a model-free method that allows us to forgo a priori assumptions on the shape and delay of the hemodynamic response function (HRF) (Di & Biswal, 2022; Hejnar et al., 2007), thus offering a more unbiased measure of the breath-hold CVR amplitude response. The PC1 loading values were compared between the AB and SCI groups for each brain region and each hemoglobin species.

2.6. Breath-hold CVR Response Delays

In addition to obtaining a measure of CVR amplitude, we investigated differences in the delays of the breath-hold CVR response by performing cross-correlation analysis. Cross-correlation analysis was applied between the preprocessed fNIRS data and the task design convolved with the canonical HRF model (Friston et al., 1998). The fNIRS time-series data were shifted by a delay range of 0 s to 20 s in increments of 0.04 s (sampling frequency of 25 Hz) with respect to the convolved task design data. This delay range was chosen based on the physiological range of breath-hold delays that has been reported in the literature (K. Chen et al., 2021; Gong et al., 2023; Moia et al., 2020; Stickland et al., 2022; Tong et al., 2011). The “optimal” delay was chosen based on the delay or lag value that yielded the highest correlation value when shifting the time-series data. The optimal delay was calculated for each subject, each channel, and each hemoglobin species. These delay values were then compared between the AB and SCI groups across all brain regions.

2.7. CVR and Injury Characteristics

After comparing the CVR amplitude and delay values between the AB and SCI groups, we further investigated whether injury characteristics in the SCI group were associated with CVR characteristics. To reduce the number of comparisons, injury characteristics were correlated with PC1 loading and delay values for any regions that showed significant differences between the AB and SCI group. For the SCI group, we investigated the duration since injury (in years) and the level of spinal injury score, calculated as the number of spinal segments affected that are below the level of injury (Table 1). The mean of the CVR amplitude and delay values were taken across all channels within the particular ROI to obtain a single ROI-specific CVR metric for each subject. The mean value was taken since not all subjects had data from all channels within a particular ROI, due to the channel pruning procedure implemented.

2.8. Statistical Analysis

We compared CVR amplitude and delays between the SCI and AB groups using mixed-effects modeling. The linear mixed-effects models were implemented in R using the lme4 package (Bates et al., 2015) and were formatted as follows: Linear mixed effects models were applied to all HbO, HbR, and HbT concentration data separately. The “Group” variable consisted of either the SCI or AB category and the “ROI” variable consisted of the 7 ROI across the fNIRS cap: the LDLPFC, RDLPFC, MPFC, LMOTOR, RMOTOR, LIPL, and RIPL regions. The “Group” and “ROI” variables were treated as fixed effects while the “Subject” variable was treated as a random effect. Not all subjects had data for all channels within a particular ROI due to the channel pruning procedure. Therefore, by implementing a mixed-effects modeling approach, we were able to account for the differences in the number of channels with sufficient data quality from each subject. All linear mixed effects models were fit by using the restricted maximum likelihood algorithm and t-tests were performed using Satterthwaite’s method (Kuznetsova et al., 2017). For analysis of associations between CVR amplitude and delays with the injury characteristics in the SCI group, Pearson’s r correlation was used for each hemoglobin species concentration.

3.1. CVR Amplitude in SCI vs. AB

Across all brain regions, there were no significant differences in CVR amplitude, as measured by PC1 Loading scores, between the SCI and AB groups (Figure 4A-C). Separate linear mixed effects models were generated for each hemoglobin concentration species and no significant differences in CVR amplitude were observed between SCI and AB groups across HbO, HbR, and HbT (Table 2). For HbT, although the AB group showed higher CVR amplitude in the right motor region, the effect was weak (fixed-effect estimate = –0.187, t = –2.214, p = 0.028). Overall, the PC1 Loading values hovered around 0.1 for both SCI and AB groups and showed similar values across hemoglobin concentration species (Figure 4).

• Download figure
• Open in new tab

Figure 4. PC1 Loadings shown for both AB and SCI groups, for all hemoglobin species:

A) oxy-hemoglobin (HbO), B) deoxy-hemoglobin (HbR), and C) total hemoglobin (HbT), across all brain regions of interest: left dorsolateral prefrontal cortex (LDLPFC), left inferior parietal lobe (LIPL), left motor cortex (LMOTOR), medial prefrontal cortex (MPFC), right dorsolateral prefrontal cortex (RDLPFC), right inferior parietal lobe (RIPL), and right motor cortex (RMOTOR). The means of the PC1 Loadings are indicated by black circles and the distribution of the AB and SCI groups are shown, denoted by green and orange, respectively.

3.2. CVR Delays in SCI vs. AB

We observed significantly longer CVR delays in the SCI group compared to the AB group in the right inferior parietal lobe, for HbR concentration (linear mixed-effects model, fixed-effects estimate = 6.565, Satterthwaite’s t-test, t = 2.663, p = 0.008) (Table 3). The average HbR CVR delay in the SCI group in the right inferior parietal lobe was 14.21 s (sd: 6.60 s), and for the AB group, the average delay in the right inferior parietal lobe was 7.08 s (sd: 7.39 s) (Figure 5). For HbO concentration, we observed slightly higher CVR delay values in the left inferior parietal lobe for the AB group compared to the SCI group (linear mixed-effects model, fixed-effects estimate = –6.027, Satterthwaite’s t-test, t = –2.518, p = 0.012) (Table 3), with the mean delay values in the AB group being 9.53 s (sd = 7.42 s) and the mean delay in the SCI group being 8.20 s (sd = 8.80 s). Similarly, for HbO concentration in the right inferior parietal lobe, we observed slightly higher CVR delay values in the AB group compared to that of the SCI group (linear mixed-effects model, fixed-effects estimate = –6.136, Sattherthwaite’s t-test, t = – 2.332, p = 0.020) (Table 3). No significant differences in delays were observed between the SCI and AB groups for HbT concentration across all ROIs (Satterthwaite’s t-tests, p > 0.05). The CVR delays varied across the different ROI investigated; however, since we were primarily interested in differences in CVR delays between the AB and SCI groups, statistical results from the linear mixed-effects models were only reported for group differences in each region (Table 3).

• Download figure
• Open in new tab

Figure 5. Breath-hold CVR delay values between AB and SCI groups for

A) oxy-hemoglobin (HbO), B) deoxy-hemoglobin (HbR), and C) total hemoglobin (HbT), across the left dorsolateral prefrontal cortex (LDLPFC), left inferior parietal lobe (LIPL), left motor cortex (LMOTOR), medial prefrontal cortex (MPFC), right dorsolateral prefrontal cortex (RDLPFC), right inferior parietal lobe (RIPL), and right motor cortex (RMOTOR). The black circles indicate the mean of the CVR delays and the distributions of the delays are shown for the AB and SCI groups, denoted by green and orange, respectively.

3.3. CVR Association with Duration Since Injury and Injury Level

Upon further investigation of the right inferior parietal region’s CVR delays, we found a significant association between the CVR delay and duration since injury for HbT concentration values in the right inferior parietal lobe. As the duration since injury increased, the CVR delay values decreased (Pearson’s r correlation, r = –0.59, p = 0.04) (Figure 6A). For HbO and HbR concentration values, no significant associations were found between the mean CVR delays in the right inferior parietal lobe and the duration since injury (Pearson’s r correlation, HbO: r = – 0.09, p = 0.82, HbR: r = 0.28, p = 0.41). Similarly, no significant associations were found between the CVR delays in the right inferior parietal lobe and the injury level score for either HbO, HbR, or HbT (Pearson’s r correlation, p > 0.05) (Figure 6B). For subjects without sufficient data in any of the channels in the right inferior parietal lobe (channels 13, 14, 17, and 18), the data were not included in the correlational analyses.

• Download figure
• Open in new tab

Figure 6. Association between the mean CVR delays in the right inferior parietal lobe and the

A) duration since injury (years) and B) spinal injury level score in the SCI group, for oxygenated hemoglobin (HbO), deoxygenated hemoglobin (HbR), and total hemoglobin (HbT) concentrations. Pearson’s correlation was performed and a least-squares line is shown in orange with the standard error of the mean in gray.

3.4. Variability in CVR Across the Brain

The SCI group showed significantly higher variability of CVR delays for HbO concentration compared to the AB group, as calculated by the standard deviation of CVR values across the brain (across all channels) (two-sample independent t-test, p = 0.0168) (Figure 7B). For CVR amplitude, there were no significant differences in variability across the brain between the AB and SCI groups across HbO, HbR, and HbT concentrations (independent-samples t-test, p > 0.05) (Figure 7A). For CVR amplitude, the whole-brain variability ranged from 0.15 to 0.25 for both AB and SCI groups, showing minimal differences in variability.

• Download figure
• Open in new tab

Figure 7. Variability in breath-hold CVR

A) amplitude (PC1 Loading) and B) delays across the brain between the AB and SCI groups, shown in green and orange, respectively. For HbO data, we found significant differences between the AB and SCI groups in the variability of delays across the brain (across channels) (two-sample t-test, p = 0.0168). Error bars indicate ±1 standard error of the mean.

4.1. No differences in CVR amplitude between the SCI and AB groups

When comparing CVR amplitude between SCI and AB groups, we observed no significant differences between the two groups. This is consistent with Weber and colleagues’ results that used fMRI to quantify CVR, in which no significant differences were observed in the standalone and steady-state CVR results between SCI and AB groups (Weber et al., 2022). They argue that unchanged full CVR may be masked by increased CVR in individuals with thoracic injuries and decreased CVR in those with cervical injuries. Furthermore, Wilson and colleagues reported that CVR is maintained in individuals with tetraplegia, whereas dynamic cerebral autoregulation may be altered, suggesting adaptive mechanisms in those with chronic tetraplegia (Wilson et al., 2010). Although our CVR metric was quantified using a hypercapnic breath-hold stimulus, the previous studies investigating CVR in individuals with SCI used CO2 gas inhalation and may be limited by the small sample sizes (7 and 6 SCI participants, respectively) (Weber et al., 2022; Wilson et al., 2010). Despite these differences, our results on the time-independent CVR metric are consistent with their conclusions. An advantage of our method to quantify CVR amplitude by using PCA is that there is no biased assumption on the specific HRF shape associated with the breath-hold CVR response, thus providing a time-independent and HRF-shape independent metric of CVR.

4.2. Individuals with SCI have significantly longer CVR delays than AB individuals

We found significantly longer CVR delays in the SCI group compared to the AB group in the right inferior parietal lobe. Longer delays in CVR is associated with poorer vascular health (K. B. Miller et al., 2019; Sobczyk et al., 2021; Stringer et al., 2021), suggesting reduced vascular health in the right inferior parietal lobe in individuals with SCI. The right inferior parietal lobe is involved in visuospatial attention and sensorimotor integration, in addition to executive function (Corbetta & Shulman, 2002; Numssen et al., 2021). Delays in this region may be due to loss of sensorimotor function in individuals with SCI after injury, particularly for the paraplegic cohort recruited in the present study. However, it is important to note that we did not find significant differences in the CVR amplitude despite finding differences in the CVR delays. This may suggest a compensatory mechanism in which a longer CVR delay in individuals with SCI is required for comparable CVR amplitudes to that of the AB group. Furthermore, a significant delay in CVR in the SCI group was observed for HbR only, which may reflect more sensitive changes to a breath-hold task in HbR rather than to HbO concentrations in the blood. CVR quantified using blood-oxygen-level-dependent (BOLD) fMRI techniques cannot fully capture differences in CVR dynamics between the different hemoglobin species. This is due to HbR being paramagnetic, causing dephasing and shortening of the transverse relaxation time (T₂*) weight, which decreases the fMRI BOLD signal; therefore, a decrease in HbR subsequently increases the strength of the signal (Mcintyre et al., 2003). HbO on the other hand, is not paramagnetic and has very little contribution to relaxation rates. Therefore, using fNIRS, we can better distinguish between CVR measured by dynamic changes not only in HbR, but also in that of HbO and HbT. This shows a unique characteristic of quantifying CVR using fNIRS as compared to that of fMRI, which relies primarily on HbR changes. Future studies may want to further investigate the different sensitivities of HbO, HbR, and HbT in quantifying CVR using fNIRS.

4.3. CVR delays in the right inferior parietal lobe is associated with the duration since injury

We further investigated the CVR delays in the right inferior parietal lobe in individuals with SCI and found a significant association between CVR delays with the duration since injury, for HbT concentration. As the duration since injury increased for the SCI group, the CVR delays decreased. This suggests that over time, the cerebrovascular health of individuals with SCI may improve, as the cerebrovasculature remodels and motor region connections reorganize after injury (Kokotilo et al., 2009). These improvements in CVR delays in the right inferior parietal lobe may also be related to the frequency of rehabilitation or exercise programs that individuals may be enrolled in. Although it is not clear the extent to which individuals with SCI in this cohort actively participate in physical rehabilitation training, generally, exercise improves cerebrovascular health, as blood perfusion is increased in the brain (Bliss et al., 2021; Davenport et al., 2012).

4.4. Variability in CVR across the brain

We found significantly higher variability in the delay of CVR across the brain in individuals with SCI compared to AB individuals, for HbO concentration. This may be explained by the differences in injury characteristics in the SCI group, whereas those in the AB group were more likely to be more homogeneous. As injury characteristics such as the completeness of injury, level of injury, duration since injury, and cause of injury differ across individuals in the SCI group, it is consistent that the CVR characteristics across the brain would also differ more in the SCI group. In healthy aging subjects, CVR has been shown to differ across brain regions, with CVR in particular regions such as the hippocampus being associated with memory scores (Catchlove et al., 2018). The heterogeneous nature of CVR has been reported in previous studies (D. Y. Chen, Di, & Biswal, 2024; D. Y. Chen, Di, Yu, et al., 2024; Kastrup et al., 1999; Moia et al., 2020; van Niftrik et al., 2023); however, it is unclear whether a measure of the heterogeneity of CVR may reveal vascular pathologies in the brain, particularly in patient populations.

4.5. Limitations

The present study contains a lack of physiological measures such as end-tidal CO2 measures and respiratory belt data. However, it is important to note that the use of a respiratory belt on patient populations such as those with SCI may be too constricting and may not be comfortable to wear while performing a breath-holding task. It is important that participants were as comfortable as possible when performing the breath-holding tasks, thus although extra physiological sensors would be helpful, it may cause other unwanted side effects from the patients. Future studies may also want to look at the differences in CVR in individuals with thoracic vs. cervical injuries and those who have tetraplegia vs. paraplegia, since autonomic dysfunction can vary greatly based on the severity of injury. Furthermore, future studies may also want to look into using breath-hold CVR as a regressor in fNIRS studies to account for possible differences between patient populations and controls that may be due primarily to vascular differences rather than neuronal differences.