A novel approach for assessing hypoperfusion in stroke using spatial independent component analysis of resting-state fMRI.

Individualized treatment of acute stroke depends on the timely detection of ischemia and potentially salvageable tissue in the brain. Using functional MRI (fMRI), it is possible to characterize cerebral blood flow from blood-oxygen-level-dependent (BOLD) signals without the administration of exogenous contrast agents. In this study, we applied spatial independent component analysis to resting-state fMRI data of 37 stroke patients scanned within 24 hr of symptom onset, 17 of whom received follow-up scans the next day. Our analysis revealed "Hypoperfusion spatially-Independent Components" (HICs) whose spatial patterns of BOLD signal resembled regions of delayed perfusion depicted by dynamic susceptibility contrast MRI. These HICs were detected even in the presence of excessive patient motion, and disappeared following successful tissue reperfusion. The unique spatial and temporal features of HICs allowed them to be distinguished with high accuracy from other components in a user-independent manner (area under the curve = 0.93, balanced accuracy = 0.90, sensitivity = 1.00, and specificity = 0.85). Our study therefore presents a new, noninvasive method for assessing blood flow in acute stroke that minimizes interpretative subjectivity and is robust to severe patient motion.

INTRODUCTION

Imaging blood flow in acute ischemic stroke is important for directing treatment decisions. In clinical practice, magnetic resonance imaging (MRI) of blood flow is typically performed using dynamic susceptibility contrast MRI (DSC‐MRI). This method is well‐established and correlates well with the in‐vivo perfusion imaging gold‐standard, 15O‐water positron emission tomography (PET; Zaro‐Weber et al., 2016; Zaro‐Weber, Moeller‐Hartmann, Heiss, & Sobesky, 2010). However, it requires the use of intravenously administered contrast agents, which are contraindicated in patients with kidney impairment (Khawaja et al., 2015) and can accumulate in tissues with repeated use (Gulani, Calamante, Shellock, Kanal, & Reeder, 2017).

Blood‐oxygenation‐level‐dependent (BOLD) functional MRI (fMRI) has been proposed as an alternative to DSC‐MRI. At rest, BOLD signal oscillations reflect a combination of different signal sources, including the hemodynamic response to spontaneous neural activity (neurovascular coupling), peripheral physiology (respiratory and cardiac function), fluctuations in blood carbon dioxide (Golestani, Chang, Kwinta, Khatamian, & Jean Chen, 2015; Wise, Ide, Poulin, & Tracey, 2004), and head motion (Liu, 2016). Disturbed cerebral blood flow is associated with changes in BOLD signal fluctuations in multiple ways. One of the best studied perfusion‐related changes in the BOLD signal is BOLD delay (Amemiya, Kunimatsu, Saito, & Ohtomo, 2014; Q. Chen et al., 2018; Christen et al., 2015; Khalil et al., 2017, 2020; Lv et al., 2013, 2018; Ni et al., 2017; Siegel, Snyder, Ramsey, Shulman, & Corbetta, 2016; Wu, Dehkharghani, Nahab, Allen, & Qiu, 2017), which are relative phase lags in low frequency oscillations (LFOs) of the BOLD signal originating from outside the brain (systemic LFOs) in underperfused brain regions (Tong, Hocke, & Frederick, 2019). Disturbed perfusion is known to affect other BOLD signal features, including a shift to lower LFO frequencies (Liu et al., 2007; Tsai et al., 2014; Yao et al., 2012) and a reduction in LFO amplitude (Wang et al., 2008). These findings suggest that regions of disturbed perfusion may have a unique set of BOLD signal characteristics that distinguish them from normally perfused tissue. Using this BOLD “signature” to rapidly and reliably identify abnormally perfused tissue would be of great diagnostic and prognostic utility.

Spatial independent component analysis (spatial ICA) is a data‐driven method for decomposing complex signals into statistically independent subcomponents (Beckmann & Smith, 2004). In fMRI, it has been used to remove nonneuronal signal contributions such as head motion, large vessel pulsations, and scanner instabilities (Thomas, Harshman, & Menon, 2002). Additionally, using spatial ICA, researchers can isolate functionally distinct resting‐state networks (RSNs) such as visual cortical areas and the sensory‐motor cortex (Beckmann, DeLuca, Devlin, & Smith, 2005).

Considering spatial ICA's ability to separate statistically independent sources of variance in the BOLD signal, we investigated whether spatial ICA would be able to detect hypoperfusion‐related changes in the BOLD signal in acute stroke patients. Specifically, we aimed (a) to identify spatial ICA components corresponding to tissue hypoperfusion, (b) to characterize the temporal and spatial properties of these “Hypoperfusion spatially independent components” (HICs), and (c) to determine whether these properties allow the automated identification of HICs.

MATERIALS AND METHODS

Patients

This study analyzes data from a subset of ischemic stroke patients from 1000Plus (Hotter et al., 2009; clinicaltrials.gov identifier: NCT00715533) and Longitudinal MRI Examinations of Patients With Brain Ischemia and Blood Brain Barrier Permeability (LOBI‐BBB, clinicaltrials.gov identifier: NCT02077582), two prospective cohort imaging studies conducted at Charité Universitätsmedizin Campus Benjamin Franklin (Berlin, Germany) between January 2009 and March 2016.

The patients included in this study had a clinical and radiological diagnosis of acute supratentorial ischemic stroke, and had received both dynamic susceptibility contrast magnetic resonance imaging (DSC‐MRI) and resting‐state functional MRI (rs‐fMRI) scans within 24 hr of symptom onset as part of a standard stroke imaging protocol (Hotter et al., 2009). All patients showed perfusion deficits on their time‐to‐maximum (T max) maps at baseline. Where available, follow‐up scans from the same patients, acquired approximately 24 hr after the baseline scans, were also included in the analysis. Exclusion criteria were any contraindications to undergoing an MRI or receiving a gadolinium‐based contrast agent.

Thirty‐seven acute ischemic stroke patients who met our inclusion criteria were included in this study, fourteen of whom had follow‐up scans. A total of 51 scans were thus included in data analysis. The baseline scans of 22 of these patients were previously analyzed in a comparison between time shift analysis of rs‐fMRI data (BOLD delay) and DSC‐MRI (Khalil et al., 2017) and 12 patients (and their follow‐up scans) from the current study were previously investigated in a study of the longitudinal evolution of perfusion deficits detected using BOLD delay (Khalil et al., 2020).

Ethics statement

Patients provided written informed consent prior to participation, and all procedures were approved by the local ethics committee (EA4/026/08 for the 1000plus study and EA1/200/13 for LOBI‐BBB) and were performed according to the Declaration of Helsinki.

Imaging

An overview of the image processing pipeline implemented in this study is shown in Figure 1.

FIGURE 1

Overview of the image processing pipeline. rs‐fMRI, resting‐state functional magnetic resonance imaging; DSC‐MRI, dynamic susceptibility contrast magnetic resonance imaging; TSA, time shift analysis; T max, time‐to‐maximum; FIX, FMRIB's ICA‐based X‐noiseifier

Acquisition

All imaging was performed on a Siemens (Erlangen, Germany) Tim Trio 3 Tesla MRI scanner. The sequence parameters were as follows, rs‐fMRI; repetition time (TR) = 2,300 ms, echo time (TE) = 30 ms, flip angle (FA) = 90°, matrix = 64 × 64, voxel dimensions = 3 × 3 × 3 mm3, 1 mm slice gap, 33 slices, and 150 volumes (acquisition time = 5 min and 50 s). The rs‐fMRI scan was performed prior to the administration of the contrast agent and patients were requested to relax, lie still, and close their eyes for the duration of the scan.

DSC‐MRI was performed after injection of a bolus of 5 ml Gadovist (Gadobutrol, 1 M, Bayer Schering Pharma AG, Berlin, Germany) followed by a saline flush at a flow rate of 5 ml/s. The sequence parameters were: TR = 1,390 ms, TE = 29 ms, FA = 60, matrix = 128 × 128, voxel dimensions = 1.8 × 1.8 × 5 mm, 0.5 mm slice gap, 21 slices, and 80 volumes (acquisition time = 1 min and 58 s).

Other imaging included diffusion‐weighted imaging (DWI), time‐of‐flight MR angiography, and fluid‐attenuated inversion recovery (FLAIR) scans, as described previously (Hotter et al., 2009).

Dynamic susceptibility contrast MRI

Maps of time‐to‐maximum of the tissue residue function (T max) were generated from the DSC‐MRI data using Stroketool (version 2.8; Digital Image Solutions—H. J. Wittsack) after deconvolution of the concentration‐time curve through block‐circulant singular value decomposition (Wu et al., 2003). The arterial input function (AIF) was selected from 5 to 10 voxels in the distal branches of the middle cerebral artery on the hemisphere contralateral to the acute infarct (Ebinger et al., 2010). Artifacts in the cerebrospinal fluid (CSF) were automatically removed from the T max maps using CSF masks derived from T2‐weighted B0 images. Perfusion deficits were manually delineated on T max maps by one of the authors (A. K., with 8 years' experience in stroke imaging research).

Resting‐state functional MRI

Preprocessing

Preprocessing of the resting‐state data was performed using tools from Analysis of Functional NeuroImages (AFNI; Cox, 1996) and the Functional Magnetic Resonance Imaging of the Brain Software Library (FSL). The first four volumes of each time series were discarded for signal equilibration. After correcting for slice‐timing effects, the remaining volumes were skull‐stripped, realigned to the mean functional image of each individual, and linear and quadratic trends were removed. Images were spatially smoothed (Gaussian kernel of full‐width‐at‐half‐maximum = 6 mm), but left temporally unfiltered to include the full range of signal frequencies for spatial ICA.

ICA decomposition

Spatial ICA is a data‐driven approach for multivariate data analysis and representation. It was implemented in this study using FSL's Multivariate Exploratory Linear Optimized Decomposition into Independent Components—MELODIC (Beckmann & Smith, 2004). Spatial ICA models four‐dimensional fMRI data as a linear combination of unknown source signals (referred to as “independent components”), each described by a spatial map and a mean time course across all voxels in the spatial map. It dissects out the different source signals by assuming their mutual, statistical independence and non‐Gaussianity (Beckmann & Smith, 2004). In this study, spatial ICA was implemented using MELODIC's default Bayesian dimensionality estimation, which estimates the number of source signals within the data automatically (Beckmann & Smith, 2004). This means that the number of independent components output by MELODIC can differ between patients and scanning sessions in the study.

Manual classification of independent components

To reduce the number of components by removing nonsignal independent components, we used FMRIB's ICA‐based X‐noiseifier (FIX), a machine‐learning algorithm that automatically classifies signal and artifact components in rs‐fMRI data using 186 spatial and temporal features (Salimi‐Khorshidi et al., 2014). FIX needs to be trained on manually classified data. Therefore, for a subset of our data (20 out of 37 baseline scans; the FIX training dataset, see below), two raters (A. A. K and J.‐Y. H.) together (i.e., not independently) manually classified the independent components from MELODIC into the following classes: resting‐state networks (RSNs), likely hypoperfusion independent components (HICs), noise (including head motion and scanner noise), or “unknown” if the independent component could not be clearly identified. Manual classification of the RSNs, noise, and unknown components was carried out according to published guidelines (Kelly et al., 2010) by inspecting the thresholded spatial maps (|Z| > 2.33), the power spectrum, and the time courses. HICs were identified as unilateral components that were present in, and largely restricted to, the vascular territory in which the patient's infarct was found (guided by the DWI).

A total of 896 independent components were generated using MELODIC from the subjects in the FIX training dataset. The median number of independent components per subject was 45 (interquartile range [IQR] = 43–47; Figure S1). The mean (±SD) contribution of each component class (% out of the total number of components) for subjects in the FIX‐training dataset are as follows: noise (79.8 ± 9.4%), resting‐state networks (9.2 ± 5.9%), and likely hypoperfusion independent component(s) (HIC; 2.6 ± 1.5%). “Other” components that could not be unambiguously classified into any of the other three categories formed 8.4 ± 4.8% of the training dataset. Of the 20 subjects in the FIX training dataset, 14 had a single HIC, 3 had two HICs, 1 had three HICs, and 2 had no HICs. Examples of our classifications are provided in Figure S2.

Data denoising

After manual classification of the FIX training dataset, minor modifications were made to FIX to accommodate our data: we substituted T1‐weighted magnetization‐prepared rapid gradient‐echo (MPRAGE) scans used in calculations of FIX features with B0 images (T2‐weighted spin‐echo EPI), as MPRAGEs are not routinely acquired as part of our acute stroke MR imaging protocol (Hotter et al., 2009). Tissue segmentation maps needed for these features were generated from these images using FSL's FMRIB Automated Segmentation Tool (FAST). All segmentations were visually assessed for quality.

FIX was trained on the FIX training dataset. The accuracy of FIX for separating signal from noise in the FIX training dataset was investigated using leave‐one‐out cross‐validation. A threshold of 20 was selected for the classification of independent components in our training data, as this was the highest threshold at which no HICs were incorrectly classified by FIX as noise (see Table S1). At this threshold, FIX achieved a mean true positive rate (TPR; percent of true signal independent components correctly classified) of 90.9%, and a mean true negative rate (TNR; percent of true noise independent components correctly classified) of 84.6%. The output of FIX denoising was a “shortlist” of components classified by FIX as being “non‐noise.” After training, FIX was applied to the full dataset (n = 2,337 independent components from 51 scans). FIX reduced the number of independent components to a list of likely‐signal independent components by 77.4% (total, n = 528 independent components from 51 scans, see Figure S3).

Time shift analysis

For a visual comparison between HICs and a more established measure of perfusion derived from rs‐fMRI, we performed time shift analysis (TSA) of the preprocessed rs‐fMRI data (Amemiya et al., 2014; Christen et al., 2015; Khalil et al., 2017, 2020; Lv et al., 2013). TSA maps were generated by assigning each voxel the value of the time shift (ranging between −20 and +20 s) that achieves maximum cross‐correlation between the voxel's time series and a recursively refined regressor derived from the global mean time series after univariate interpolation of the time series. A long tracking range was applied because time shift delays in acute stroke patients are very prolonged (Khalil et al., 2017, 2020; Tanrıtanır et al., 2020) TSA was performed using a set of Python tools (https://github.com/bbfrederick/rapidtide; Frederick, 2016).

Image registration

All images (including T max maps, resting‐state data, and spatial ICA maps) were registered to a custom echo planar imaging (EPI) template in Montreal Neurological Institute 152 (MNI152) standard space using a rigid body spatial transformation (3 rotations and 3 translations, implemented using FMRIB's Linear Image Registration Tool—FLIRT) for further processing. This template was derived from the EPI scans of 103 stroke patients and details of how it was generated can be found in the supporting information of Khalil et al. (2017). This template was used instead of standard templates, such as the MNI template, to account for differences in brain and ventricular system size between the study's population and younger, healthy individuals.

Feature extraction

We extracted a set of relevant spatial and temporal features from the shortlisted independent components (i.e., after data denoising by FIX) to examine how HICs may differ from other signal components. The following features were extracted:

Temporal delay to a reference signal: Temporally delayed BOLD signals, calculated using time shift analysis (Lv et al., 2013), have been found to correspond to hypoperfused brain regions. We thus explored whether the temporal delay between the time courses of each independent component and a reference time course was different between HICs and independent components representing normally perfused tissue. We used two different reference time courses: the whole brain (global) signal (Amemiya et al., 2014; Lv et al., 2013) and the signal from within the major venous sinuses (Aso, Jiang, Urayama, & Fukuyama, 2017; Christen et al., 2015; Khalil et al., 2017), the latter of which was automatically extracted using a venous sinus template in MNI space [details described in the supporting information of Khalil et al., 2017]. Positive delay values indicate the component's BOLD signal time course followed the reference time course, while negative delay values indicate the component's BOLD signal time course preceded the reference time course.

Power in frequency bands: Ischemia shifts BOLD signal oscillations toward lower frequencies (Liu et al., 2007; Tsai et al., 2014; Yao et al., 2012). For each independent component, we calculated the percentage of total spectral power in six different frequency ranges (0–0.01 Hz, 0.01–0.025 Hz, 0.025–0.05 Hz, 0.05–0.1 Hz, 0.1–0.15 Hz, and 0.15–0.20 Hz) after fast Fourier transform of the time courses (Salimi‐Khorshidi et al., 2014).

Percent restriction to a single vascular territory: Ischemic stroke‐related hypoperfusion is most commonly restricted to a vascular territory supplied by a single major artery. This makes the percent restriction of an independent component's spatial map to a single vascular territory a potential feature of interest. We thus calculated the percentage of each independent component's volume that exists within a single vascular territory. For this analysis, only independent components with >50% of their volume within a single vascular territory were considered. The vascular territories were defined according to a custom vascular territory atlas in MNI template space, described in further detail in the supporting information of Khalil et al. (2017).

Mean T max delay: To compare the results of the spatial ICA with the delay in perfusion revealed by DSC‐MRI in the area covered by each independent component, we calculated the mean T max delay of all voxels in the T max map overlapped by the independent component's spatial map.

For each patient, min–max normalization was applied to the values of the features extracted from each of the independent components using the following equation: where, z is the normalized feature value, x is the original feature value, i refers to the independent component in question, and S refers to the feature values of the entire set of independent components from a single patient.

Automated identification of HICs

Both model training and testing were performed on the data after excluding follow‐up scans.

Model training

We used a generalized linear model (GLM) on the extracted features to estimate the probability that a given independent component is a HIC. To mitigate multicollinearity effects we applied elastic net regularization to our GLM (Zou & Hastie, 2005).

Disregarding components classified as noise, approximately 10% of the components in each dataset were HICs. To address this problem of extreme class imbalance (Krawczyk, 2016), which can bias GLM results, we first sub‐sampled the FIX training dataset to contain 50% HICs (number of components = 23), and 50% non‐HICs (number of components = 23).

For the model training, components were classified as HICs if their model probabilities were above a probability threshold determined by Youden's index, which gives equal weight to sensitivity and specificity. The shrinkage parameter lambda was estimated using leave‐one‐out cross‐validation on the training dataset.

Model testing

The test dataset contained all subjects who were not part of the FIX training dataset (number of HICs = 22, number of non‐HICs = 151). From this, we randomly subsampled 5 HICs and a set of non‐HICs according to the mean ratio between HICs and non‐HICs in individual MELODIC outputs in the FIX shortlist (about 1:10). Note that this allowed the testing of model performance to be performed on a dataset with a class imbalance similar to that of the original dataset. This random selection was repeated 50 times and the model's performance is reported as the median of these iterations (see below for metrics of model performance).

Statistical analysis

Radar plots were created using the “radarchart” function from the fmsb R package (Nakazawa, 2018) to visualize feature values. Raincloud plots, which combine dot plots, box plots, and violin plots, are used to visualize the distribution of continuous variables (Allen, Poggiali, Whitaker, Marshall, & Kievit, 2019).

The elastic net regularized GLM used in this study was implemented using the “cv.glmnet” function from the glmnet R package (Friedman, Hastie, & Tibshirani, 2010).

Model performance is reported using the following metrics:

Area‐under‐the‐curve of the receiver operating ROC curve (AUC): The probability that the model will assign a higher probability of being a HIC to a randomly chosen HIC than to a randomly chosen non‐HIC. ROC curves were visualized and their AUCs were calculated using the pROC R package (Robin et al., 2011).

Sensitivity: The true positive rate, or the proportion of HICs correctly classified by the model as being HICs.

Specificity: The true negative rate, or the proportion of non‐HICs correctly classified by the model as being non‐HICs.

Balanced accuracy: (Sensitivity + specificity)/2. Balanced accuracy was chosen instead of regular accuracy due to the unbalanced nature of our test dataset (Brodersen, Ong, Stephan, & Buhmann, 2010).

Cohen's kappa: The agreement between the model's classification and the “true”/”reference” classification, beyond that expected due to random chance. This was calculated using the irr R package (Gamer, Lemon, Fellows, & Sing, 2019).

The GLM coefficients were converted into Odds ratios (Odds ratio = e coefficient). Since elastic net, like all penalized regression models, provides biased estimators, no meaningful standard errors can be calculated (Kyung, Gill, Ghosh, & Casella, 2010) and therefore P‐values are not reported in this analysis. Instead, the model produces penalized coefficients through bootstrapped cross‐validation. Variables whose contribution to the outcome is negligible have penalized coefficients of zero (i.e., Odds ratios of 1).

Data/code availability

The data and code for the training and testing of the model can be found at https://github.com/ahmedaak/spatial_ICA_stroke.

RESULTS

The study sample consisted of 17 women and 20 men (mean age = 71 years, SD = 14 years). Their median National Institutes of Health Stroke Scale (NIHSS) was 7 (IQR 3–14) and the median time from symptom onset to MRI was 7 hr (IQR 2–16 hr). Intravenous thrombolysis was administered to 19 (51%) patients. Fourteen patients received follow‐up scans on the second day following stroke onset, resulting in a total of 51 scans. Tissue reperfusion was successful in 3 of the 14 patients with follow‐up scans, as demonstrated by the complete resolution of time‐to‐maximum (T max) lesions observed on their baseline scans.

Head motion was quantified using framewise displacement (FD), which is the sum of the absolute value of individual subjects' six translational and rotational realignment parameters (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012). Twelve out of the fifty one rs‐fMRI scans exhibited severe head motion, defined according to previous studies on the use of rs‐fMRI for assessing cerebral hemodynamics as a mean FD of >0.4 mm across the scan or a maximum FD of >3 mm (Khalil et al., 2017, 2020). Median FD (mean across all scan volumes) in the sample was 0.23 mm (IQR: 0.15 mm) for rs‐fMRI, and 0.29 mm (IQR: 0.19–0.47 mm) for DSC‐MRI. Median FD (maximum across all scan volumes) were 1.28 mm (IQR: 0.61–2.86 mm) and 0.82 (IQR: 0.49–1.89 mm) for rs‐fMRI and DSC‐MRI scans, respectively. The distribution of mean and maximum FD for the rs‐fMRI scans is shown in Figure S4.

Hypoperfusion independent components reflect tissue hypoperfusion

Of the 37 baseline scans, 34 showed at least one HIC (6 scans showed 2 HICs and 2 scans showed 3 HICs). Figure 2 shows the DWI, T max, and HIC of three example patients. Maps for the other 31 patients can be found in Figure S5. Together, these figures reveal the striking spatial similarities between HICs and perfusion deficits on T max maps. The spatial overlap between HICs and manually delineated T max perfusion deficits was quantified using the Dice similarity coefficient (Dice, 1945). Figure S6 shows that the HIC was the independent component that showed the highest spatial overlap with the T max perfusion deficit in 40 out of the 45 datasets in which a HIC was identified.

FIGURE 2

Spatial independent component analysis (spatial ICA) of resting‐state fMRI detects post‐stroke perfusion deficits in the form of hypoperfusion independent components (HICs). The figure shows (from left to right) diffusion‐weighted imaging (visualization of infarcted tissue), hypoperfusion independent components (HICs), and time‐to‐maximum of the residual curve from dynamic susceptibility contrast MRI (T max; reflecting delayed perfusion), from three representative acute stroke patients with a right‐sided posterior cerebral artery infarct (patient A), right‐sided middle cerebral artery infarct (patient B), and left‐sided middle cerebral artery infarct (patient C). A corresponding perfusion deficit is seen in the affected vascular territory of each patient on the T max maps. In each case, spatial ICA reveals a HIC with striking spatial resemblance to the perfusion deficit visible on the T max maps obtained using a contrast agent

In addition, 10 out of the 14 follow‐up scans showed HICs. Figure 3 illustrates the relationship between the spatial ICA results and vessel status (recanalization or persistent vessel occlusion) in two patients with follow‐up scans. In patient A, we observe the resolution of HICs following vessel recanalization and tissue reperfusion, while in patient B, HICs persist with continued vessel occlusion and delayed tissue perfusion. This finding suggests that HICs dynamically respond to changes in vessel status and perfusion dynamics over time.

FIGURE 3

Perfusion assessed using spatial independent component analysis (spatial ICA) evolves with vessel status. Baseline images were obtained within 24 hr of symptom onset, while follow‐up images were acquired the next day. Patient A had a left‐sided acute posterior cerebral artery (PCA) infarct due to occlusion of the P2 segment of the PCA, seen on time‐of‐flight magnetic resonance angiography (TOF‐MRA) and diffusion weighted imaging (DWI). A perfusion deficit as seen on a contrast‐enhanced T max map and a corresponding hypoperfusion independent component (HIC) were observed among 13 other independent components (after automatic denoising) at baseline. Following recanalization of the PCA, there is a normalization of the perfusion deficit in the PCA vascular territory on the contrast‐enhanced T max map. Spatial ICA on follow‐up imaging generated 20 independent components after automatic denoising, mostly corresponding to resting‐state networks, and no HICs. Patient B had a right‐sided acute PCA infarct due to occlusion of the P1 segment. A perfusion deficit and corresponding HIC are seen at baseline. Spatial ICA also produced 10 other independent components that mostly corresponded to resting‐state networks (RSNs). Unlike in patient A, patient B's PCA remained occluded at follow‐up. The perfusion deficit in the posterior cerebral artery territory persisted on the T max map, and a HIC was observed along with RSN independent components. This demonstrates that HICs follow longitudinal perfusion dynamics corresponding to clinically relevant events such as changes in vessel status

HICs were identified in 10 out of 12 subjects with high head motion during their rs‐fMRI scans, as shown in an example patient in Figure 4. The images and head motion traces of all 12 subjects with high head motion can be interactively viewed here: https://doi.org/10.6084/m9.figshare.13676779. Overall, there were 4 false negative cases, in which no HIC was visible despite the presence of a perfusion deficit on T max maps (these cases can be interactively viewed here https://doi.org/10.6084/m9.figshare.13664384). In two of the false negative cases, extremely severe motion was present in the rs‐fMRI data, exceeding 20 mm maximum framewise displacement (see Figure S4). In the other two false negative cases, motion during the rs‐fMRI scan was much lower (mean FD = 0.25 and 0.14 mm, max FD = 0.56 and 0.44 mm).

FIGURE 4

Hypoperfusion independent components (HICs) are detected using spatial ICA despite high head motion. In this patient with a right‐sided middle cerebral artery infarct (red circles on the diffusion‐weighted image), the top row demonstrates the framewise displacement of the patient's head over the duration of the resting‐state functional MRI scan (maximum motion = 3.2 mm, mean motion = 0.35 mm). Despite the high motion, a HIC was detected that corresponded to the patient's perfusion deficit as detected by contrast‐enhanced T max maps on dynamic susceptibility contrast MRI. No perfusion deficit is seen on the time shift analysis map, derived from the same scan as the HIC. Of note, this patient has two distinct perfusion deficits within the same vascular territory (right middle cerebral artery)—one frontal and one temporoparietal—both of which are reflected in the HIC

HICs identified multiple perfusion deficits located in different vascular territories within the same patient (Figure 5). The fact that multiple, spatially remote areas of delayed perfusion within an individual patient are captured in a single independent component suggests that these areas share a BOLD signal signature that sets them apart from other components.

FIGURE 5

Hypoperfusion independent components (HICs) detect multiple perfusion deficits affecting different vascular territories in the same patient. This patient has multiple infarcts (with corresponding perfusion deficits on the contrast‐enhanced T max maps—bottom row)—in the left posterior cerebral artery (PCA) territory (right) and in the right posterior inferior cerebellar artery (PICA) territory (left). Spatial independent component analysis reveals a single HIC (middle row) that corresponds to both perfusion deficits in this patient. This implies that these areas of delayed perfusion, although spatially distanced from each other, share a common BOLD signal signature. Note that, in the left column, the DWI slice showing the PCA infarct does not correspond to the same slice on the T max and HIC images, as the infarct is not visible on the DWI at the level shown on the T max and HIC images

HICs show unique characteristics

Table 1 shows the feature values (mean across patients) extracted from the independent components. The means of the normalized feature values across patients are shown on the radar plot in Figure 6. HICs have a unique feature signature compared with non‐HICs, characterized by delayed low frequency oscillations compared with the whole‐brain and venous sinus references, greater restriction to a single vascular territory, greater overlap with areas of delayed perfusion (as assessed using T max), and more power in the lowest frequency bands (0–0.025 Hz).

TABLE 1

Mean across patients and 95% confidence interval of the feature values (non‐normalized) for HICs and non‐HICs

Feature	HICs	Non‐HICs
Delay to whole‐brain signal (s)	7.2 (5.4 to 9.1)	−0.5 (−0.8 to ‐0.1)
Delay to venous sinus signal (s)	5.9 (4.1 to 7.6)	−2.2 (−2.6 to ‐1.9)
T max delay (s)	2.15 (1.8 to 2.5)	1.3 (1.2 to 1.4)
Vascular territory occupancy (%)	51.6 (44.3 to 58.9)	6.2 (4.6 to 7.8)
Power in frequency bands (%)
0–0.01 Hz	21.7 (19.1 to 24.3)	17.9 (17.3 to 18.5)
0.01–0.025 Hz	23.9 (22.1 to 25.7)	21.4 (20.9 to 22.0)
0.025–0.05 Hz	22.6 (21.2 to 24.0)	25.8 (25.4 to 26.3)
0.05–0.1 Hz	15.1 (13.6 to 16.7)	17.4 (17.0 to 17.9)
0.1–0.15 Hz	12.5 (11.2 to 13.8)	13.3 (12.9 to 13.8)
0.15–0.2 Hz	4.2 (3.6 to 4.8)	4.0 (3.8 to 4.2)

FIGURE 6

Hypoperfusion independent components (HICs) are characterized by distinct features that distinguish them from non‐HICs independent components. This radar plot shows the feature values for HICs (n = 57) and non‐HICs (n = 471). Solid lines represent the mean feature value and dashed lines represent the 95% confidence interval around the mean. Compared with non‐HICs within the same patient, HICs have specific characteristics: higher percentage occupancy within a single vascular territory, presence in regions of higher T max delay, time courses that show a higher temporal delay to the whole‐brain and the venous sinus time courses, and higher signal power in the lowest frequency bands (0–0.01 Hz and 0.01–0.025 Hz). This combination of features represents a unique signature for HICs

Hypoperfusion independent components can be automatically identified

Elastic net regularized GLM was applied to the FIX‐shortlisted independent components (i.e., the independent components that were classified by FIX as being non‐noise) to determine the predictive power of a set of features on HIC classification. The model discriminated HICs from non‐HICs with a median balanced accuracy of 0.90, sensitivity of 1.00, and specificity of 0.85. The median Cohen's kappa for agreement between the model's classification and the reference classification was 0.51 and the median AUC was 0.93. Figure 7 shows the distribution of these performance metrics across 50 iterations of the model, using different combinations of HICs and non‐HICs.

FIGURE 7

Hypoperfusion independent components (HICs) can be automatically distinguished from other independent components using machine‐learning. Using independent test data of different combinations of HICs and non‐HICs (that were not used for model training), we ran 50 iterations of the elastic net regularized generalized linear model (GLM) to classify HICs. The overall model performance is shown in the receiver operating characteristic (ROC) curve (a). Each iteration of the model is shown in blue and the mean ROC curve is shown as a black dashed line. The ROC curve depicts the true positive and false positive rates for each classification threshold in the model. The raincloud plot in (b) shows the distribution of different metrics of model performance across the 50 model iterations. The lower and upper hinges of the box plot represent the 25th and 75th percentiles, respectively, and the horizontal bar represents the median value. Across all model iterations, the values of area‐under‐the‐ROC‐curve (AUC), balanced accuracy, sensitivity, specificity, and Cohen's kappa are depicted

A useful feature of elastic net regularized GLM is that it performs implicit variable selection, and thus enhances model interpretability by allowing the assessment of each feature's relative importance for the classification (Zou & Hastie, 2005). Table 2 shows that 5 out of the 9 extracted features were used by the model to distinguish HICs from non‐HICs (i.e., have an Odds ratio ≠ 1). These features included a larger delay in the low frequency oscillations of the HIC relative to the whole‐brain and venous sinus signals, a larger restriction to a single vascular territory by the HICs, and a higher power in the 0.025–0.05 Hz and 0.1–0.15 Hz frequency bands in the non‐HICs.

TABLE 2

Results of the elastic net regularized generalized linear model for the classification of HICs based on spatial and temporal features

Feature	Odds ratio
Delay to whole brain signal	1.14
Delay to venous sinus signal	1.19
Vascular territory occupancy	1.02
Power in frequency bands
0–0.01 Hz	1.00
0.01–0.025 Hz	1.00
0.025–0.05 Hz	0.07
0.05–0.1 Hz	1.00
0.1–0.15 Hz	0.13
0.15–0.2 Hz	1.00

DISCUSSION

In this study, we present a novel approach for identifying areas of disturbed perfusion in ischemic stroke without the use of contrast agents. We demonstrate that, using this method, areas of delayed perfusion can be extracted from resting‐state fMRI data in a manner that is user‐independent and robust to patient motion. The results of this study have important implications from both a clinical and methodological standpoint.

We provide several lines of evidence that certain independent components isolated from rs‐fMRI data using spatial ICA reflect regions of delayed brain tissue perfusion. Firstly, there is high visual similarity between the spatial maps of HICs and the T max maps of individual patients. This is noticeable even when patients had multiple perfusion deficits within (Figure 4) and across (Figure 5) vascular territories. Secondly, HICs were present in areas of much higher T max delay than other independent components (Figure 6). Thirdly, HICs observed at baseline disappeared at follow‐up in patients with successful recanalization and reperfusion, paralleling the reversal of patients' DSC‐MRI T max lesions (Figure 3). Finally, HICs show BOLD signal characteristics that are consistent with those previously described in hypoperfused tissue (Liu et al., 2007; Lv et al., 2013; Tsai et al., 2014; Yao et al., 2012).

These spatial and temporal features of HICs help distinguish them from resting‐state networks (Figure 6). BOLD oscillations associated with HICs were more temporally delayed in relation to the global and venous sinus signals than oscillations associated with resting‐state networks. Although temporal delays in BOLD LFOs have been shown to occur under physiological conditions (Aso et al., 2017) and differ between resting‐state networks (Tong, Hocke, Fan, Janes, & Frederick, 2015), the largest delays are seen in areas of pathologically reduced perfusion (Khalil et al., 2017, 2020). Additionally, HICs showed more signal power in very low frequency ranges (0–0.025 Hz) and less signal power in the 0.025–0.05 Hz range than resting‐state networks. This finding is consistent with previous studies showing increased signal power in the <0.01 Hz range in hypoperfused tissue (Liu et al., 2007; Tsai et al., 2014; Yao et al., 2012). Our study adds to this by finding that oscillations associated with resting‐state networks exist largely toward the higher end of the low frequency range than oscillations associated with disturbed perfusion.

One disadvantage of spatial ICA is that it can output potentially dozens of components, and browsing through them to visually identify those likely reflecting hypoperfusion can be time‐consuming, require expertise, and be subject to bias. We therefore combined the aforementioned component features with the degree of restriction of the component's spatial map to a single vascular territory and used this set of features to train an algorithm to automatically distinguish HICs from resting‐state networks. The algorithm did this with high (>95%) balanced accuracy, sensitivity and specificity, showing that HICs can be automatically extracted from the rest of the components of spatial ICA.

In this study, the spatial agreement between HICs and delayed perfusion is reflected in the fact that the areas covered by the HICs on T max maps showed by far the highest T max delay (after removal of artifactual T max delays in the CSF). In addition, in the majority of datasets, the highest spatial overlap was found between the manually delineated T max perfusion deficits and the HICs (Figure S6). Finally, the visual similarity between HICs and T max maps (which can be interactively viewed here: https://doi.org/10.6084/m9.figshare.13686931), provides further evidence for their spatial agreement.

Spatial ICA has certain advantages over existing methods for assessing perfusion using rs‐fMRI. Instead of assessing the individual characteristics of hypoperfusion‐related signals separately, spatial ICA extracts components that reflect a combination of changes in frequency, amplitude, and temporal delay. In addition, we show that spatial ICA is capable of identifying perfusion deficits even in the presence of severe patient motion, where methods based on cross‐correlation, such as time shift analysis (Lv et al., 2013), often fail (Figure 4 and https://doi.org/10.6084/m9.figshare.13676779.v1). Finally, unlike time shift analysis, spatial ICA does not require the specification of a reference signal, the choice of which can substantially affect the calculated maps (Christen et al., 2015; Khalil et al., 2017; Wu et al., 2017).

This study's findings have two major implications. From a clinical perspective, it provides a proof‐of‐concept for a new method of assessing blood flow in acute stroke patients that may be relevant for informing decision‐making regarding recanalization and reperfusion therapies. This method is safer, as it does not require the use of exogenous contrast agents, which are particularly problematic in this patient population, who are generally older and have a higher prevalence of chronic kidney disease (Sadowski et al., 2007). The acquisition of the rs‐fMRI data, upon which this method is based, is simple and widely available on clinical scanners. Although the rs‐fMRI data acquisition took substantially longer than DSC‐MRI in this study (about 6 vs. 2 min), recent evidence suggests that extracting perfusion‐related information from rs‐fMRI can be achieved with much shorter acquisitions (Tanrıtanır et al., 2020). The method is also user‐independent, thereby avoiding the subjective assessment of perfusion maps that often leads to inconsistencies between and within experts (Campbell et al., 2010) and allows it to be seamlessly integrated in routine clinical practice. Finally, the robustness of the method to head motion is a substantial practical advantage, as patients scanned during acute illnesses tend to exhibit a lot of motion, which diminishes scan quality and interpretability (Andre et al., 2015). For further clinical validation of this method, detailed comparisons to quantitative perfusion thresholds derived from DSC‐MRI and to imaging and clinical outcomes should be made in future studies.

From a methodological perspective, this study underscores the importance of accounting for disturbed perfusion as a potential source of confounding in rs‐fMRI studies, particularly in patients with cerebrovascular diseases. More consideration should be given to disentangling the effects of disturbed perfusion on the BOLD signal (e.g., due to vessel pathology) from components of the BOLD signal that reflect neuronal activity via local neurovascular coupling. So far, suggestions on how to do this have included regressing out the time delays (relative to a reference) from the rs‐fMRI data (Erdoğan, Tong, Hocke, Lindsey, & De Frederick, 2016) and temporally realigning the BOLD signal time courses according to each voxel's time delay value (Jahanian, Christen, Moseley, & Zaharchuk, 2018). Considering, however, that tissue hypoperfusion is a pathophysiological hallmark of stroke, the decision to remove the effect of disturbed perfusion from the BOLD signal will depend on the exact research question.

While the spatial distribution of physiological vascular processes overlap with, and are often indistinguishable from, the spatial distribution of resting‐state networks (Bright, Whittaker, Driver, & Murphy, 2020; Chen et al., 2020; Tong et al., 2015), we show that the influence of disturbed perfusion on the BOLD signal in stroke patients can be readily disentangled from other components of the BOLD signal using spatial ICA. There are two potential explanations for this. The first is that, in stroke, hypoperfusion is spatially restricted to either a vascular territory or part of a vascular territory. On the other hand, physiological vascular processes are spatially distributed in a manner similar to resting‐state networks (Bright et al., 2020; Chen et al., 2020; Tong et al., 2015), which may make their separation using spatial ICA less likely. The second possible explanation has to do with the large difference in the temporal BOLD signal characteristics between hypoperfused and normally perfused tissue, which may also facilitate their separation by spatial ICA. In stroke, disturbed perfusion tends to be severe compared with the physiological delays in perfusion across different brain regions, and this leads to relatively large changes in the temporal characteristics of the BOLD signal in hypoperfused regions (Khalil et al., 2017, 2020; Liu et al., 2007; Tsai et al., 2014; Yao et al., 2012).

As the first study to describe this method for assessing perfusion, the study has some limitations. The study sample is relatively small, owing to the fact that it includes an established, albeit relatively invasive, reference standard for assessing perfusion (DSC‐MRI) as a comparison. Because of the novelty of the method, we chose to have the two raters perform the ratings together and therefore could not quantify interrater agreement on the identification of HICs in this study. Future studies should test the algorithm we developed on larger cohorts with different MR sequence parameters and more heterogeneous patient cohorts, assess practical aspects of the method, such as the required computing time and power, and investigate the minimum scan length required for the method to deliver reliable results as has recently been done for time shift analysis (Tanrıtanır et al., 2020). Larger cohorts would also allow the investigation of the clinical significance of perfusion assessed using this method, in terms of how it relates to clinical and imaging outcomes and how it potentially influences clinical decision‐making in acute stroke patients (Fisher & Albers, 2013). Only two of the four cases in our cohort where spatial ICA could not identify the perfusion deficit could be explained by severe head motion. Therefore, larger studies should investigate in more detail the causes of such false negative cases. Finally, the nature of ICA means that, depending on the properties of the algorithm used, individual signal sources (such as a resting state network or hypoperfused tissue) can be spread across multiple components (Esposito et al., 2002). Therefore, it is important to note that the term “hypoperfusion independent component” is an operationalization and that other independent components might also partially reflect hypoperfused tissue.

In summary, spatial independent component analysis is a novel approach for identifying hypoperfused tissue in ischemic stroke. It does not require the use of exogenous contrast agents, its data can be analyzed without user input, and its results are robust to patient motion. It therefore presents a convenient and promising new alternative to existing perfusion imaging methods in acute stroke.

CONFLICT OF INTEREST

J‐Y.H., A.A.K., K.V., and J.B.F. are co‐inventors of a method for automatically delineating perfusion lesions on perfusion MRI data (European Patent application 17179320.01‐1906), distinct from the method described in this study, which is now in the public domain. J.B.F. reports grants from European Union Seventh Framework Program and personal fees from Bioclinica, Artemida, Cerevast, Brainomix, BMS, Merck, Eisai, Biogen, Guerbet, and Nicolab outside the submitted work. All other authors declare that they have no conflict of interest.

Figure S1 Distribution of independent component classes in the training dataset. Bar plot showing the results of the manual classification of independent components for individual patients in the training dataset (n = 20). Each bar represents one rsfMRI scan from one patient. The total number of components differs across patients because MELODIC estimates data dimensionality automatically on a per‐subject basis. RSN, resting‐state network; HIC, hypoperfusion independent component.

Figure S2 Examples of manual classifications of the output of spatial independent component analysis (spatial ICA). Each independent component was classified as either resting state network, hypoperfusion independent component, noise (subdivided into arteries, veins, head motion, or scanner), or other (when unclear). Classification was based on the spatial distribution of the component map, as well as its temporal characteristics.

Figure S3 Results of FMRIB's ICA‐based X‐noiseifier (FIX) denoising of the full dataset. The bar plot shows the number of likely‐signal independent components (pink) identified by FIX from the total number of components (gray) for each subject (x‐axis) in the full dataset of 51 scans.

Figure S4 Distribution of head motion during the rs‐fMRI scan in the study sample. Scatter plot showing the distribution of mean and maximum framewise displacement values in the study sample. The data points marked by a red cross correspond to patients in which no hypoperfusion independent component (HIC) was visible despite the presence of a perfusion lesion on T max maps (i.e., false negatives). Note that both the x and y axes are on a logarithmic scale.

Figure S5 Spatial independent component analysis (spatial ICA) of resting‐state fMRI detects post‐stroke perfusion deficits in the form of hypoperfusion independent components (HICs). The figure shows diffusion‐weighted imaging (visualization of infarcted tissue), time‐to‐maximum of the residual curve (T max; reflecting perfusion), hypoperfusion independent components (HIC), and time shift analysis (TSA) maps from all patients in the study sample (n = 37). The top part of the figure shows patients who only received baseline scans (n = 23), while the bottom part shows patients who received both baseline and follow‐up scans (n = 14). An “N” in the HIC panel represents a false‐negative; absent HIC despite the presence of a corresponding T max lesion. An “R” in the HIC panel (only in the lower part of the figure) represents an absent HIC due to vessel recanalization and tissue reperfusion on the follow‐up scans. The dashed red boxes around some patients indicate severe head motion during the rsfMRI scan (>0.4 mm mean framewise displacement or >3 mm maximum framewise displacement). Note that, for the time shift analysis maps, only “positive delays” indicating hypoperfused tissue are shown for easier visualization.

Figure S6 Spatial overlap between the independent components and the T max perfusion deficits. Each scatter plot represents a dataset from the study, excluding the three follow‐up datasets where reperfusion had occurred and therefore no perfusion deficit was visible on the T max maps. The y‐axis shows the Dice similarity coefficient—a measure of spatial overlap ranging from 0 (no overlap) to 1 (perfect overlap) between the manually delineated T max perfusion deficits and each independent component from each dataset. The x‐axis shows the index of each independent component in the dataset. Red points represent hypoperfusion independent components (HICs). Note that in three cases, no HIC was identified despite the presence of a T max perfusion deficit (false negatives). The HIC was the independent component that showed the highest spatial overlap with the T max perfusion deficit in 40 out of the 45 datasets in which a HIC was identified.

Table S1 FMRIB's ICA‐based X‐noiseifier (FIX) classification accuracies from training data (n = 20).

Click here for additional data file.