Analysis of dynamic texture and spatial spectral descriptors of dynamic contrast-enhanced brain magnetic resonance images for studying small vessel disease

Cerebral small vessel disease (SVD) comprises various pathological processes affecting small brain vessels and damaging white and grey matter. In this paper, we propose a framework comprising region of interest sampling, dynamic spectral and texture description, functional principal component analysis, and statistical analysis to study exogenous contrast agent distribution over time in various brain regions in patients with recent mild stroke and SVD features.We compared our results against current semi-quantitative surrogates of dysfunction such as signal enhancement area and slope. Biological sex, stroke lesion type and overall burden of white matter hyperintensities (WMH) were significant predictors of intensity, spectral, and texture features extracted from the ventricular region (p-value < 0.05), explaining between a fifth and a fourth of the data variance (0.20 ≤Adj.R2 ≤ 0.25). We observed that spectral feature reflected more the dysfunction compared to other descriptors since the overall WMH burden explained consistently the power spectra variability in blood vessels, cerebrospinal fluid, deep grey matter and white matter. Our preliminary results show the potential of the framework for the analysis of dynamic contrast-enhanced brain magnetic resonance imaging acquisitions in SVD since significant variation in our metrics was related to the burden of SVD features. Therefore, our proposal may increase sensitivity to detect subtle features of small vessel dysfunction. A public version of the code will be released on our research website.


Introduction
Cerebral small vessel disease (SVD) encapsulates multiple pathological processes disrupting the optimum functioning of perforating cerebral arterioles, capillaries, and some venules, resulting in grey matter (GM) and white matter (WM) damage [1][2][3]. SVD is a serious problem causing between 20% to 25% of strokes, up to 45% of dementias, and substantial cognitive, psychiatric, and physical disabilities. At a global scale, SVD may be leading to between three and four million new cases of stroke [4] and 16 million new cases of dementia per year 1 [3]. Despite being a worldwide matter and governmental priority, little is known about its cause(s) since much of SVD is clinically silent and late [1,4]. Therefore, efforts for understanding the pathophysiological mechanisms surrounding SVD and developing techniques to characterise this disease are of global impact.
Dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) is commonly used to investigate endothelial dysfunction, a pathophysiological component thought to be associated with the SVD pathogenesis [6], as it permits detecting leakage in tissue and cerebrospinal fluid (CSF) spaces thought to be caused by an impaired bloodbrain barrier (BBB) or blood-CSF barrier [7][8][9]. In this imaging modality, a series of acquisitions are taken before and after injecting a Gadolinium-based contrast agent intravenously to image its distribution through brain tissues over time. The contrast agent causes the relaxation time of water molecules to decrease in T1w. Therefore, its accumulation in the extracellular extravascular space with increased BBB impairment leads to increased signal enhancement.
In a recent work [10], and CSF cavities were quantitatively measured in pre-and post-contrast study and showed to vary with increased overall SVD and white matter hyperintensity (WMH) burden. In this work, we aim to study signal intensity fluctuation on the entire DCE-MRI sequence using established computer vision descriptors to determine whether specific variation patterns relate to the health of the patient and their suitability for acting as a surrogate measure of small vessel impairment. We consider texture and spectral descriptors. Textures, which encode local signal changes within a region of interest, have been successfully used to characterise WMH in T2 FLAIR scans [11]; study prevs post-contrast differences in small vessel disease [10]; and classify breast lesions into benign and malignant, predict chemotherapy response, and diagnose prostate cancer [12] in DCE-MRI. Power spectra, the strength of frequency components into the overall signal, have been successfully applied in dynamic susceptibility contrast MRI to characterise neurophysiological and hemodynamic patterns of Alzheimer's disease [13], detect and characterise oscillations in blood oxygen leveldependent imaging reflecting network connectivity [14], and discern between conduct disorder and healthy subjects from resting functional MRI acquisitions [15]. We hypothesise that the application of these computer vision descriptors to DCE-MRI scans can identify tissue differences in brain regions that relate to the burden of SVD features.
We propose a framework to study contrast signal-time trajectory in healthy and pathological brain regions in DCE-MRI acquisitions. The proposal comprises region of interest sampling, dynamic texture and spatial spectral description, functional principal component analysis, and statistical group comparison and linear regression. The contributions of this work are: (i) we introduce a fully functional framework to analyse DCE-MRI acquisitions based on dynamic spectral and texture descriptors and (ii) we show an application of our framework to the study of DCE-MRI signals of a cohort (n = 42) with a wide range of SVD burden.

Materials and methods
The pipeline consists of four steps, as shown in Fig. 1. First, we segmented all regions of interest for each patient in the cohort. Second, we described signals using dynamic spectral and texture descriptors. Third, we examined the resulting descriptions using a functional principal component analysis (FPCA). Fourth, we studied whether scores of the primary mode of variation were associated with any of the clinical variables. Details of each step are provided in the following sections.

Subjects and clinical variables
We used data from a prospective study of patients with recent mild stroke and SVD features (n = 42 subjects, 12 women, 19 lacunar stroke). Of 201 in the original study, we selected 42 on the basis of considering only high-quality scans (qualitative assessment of truncation and motion artefacts) and representing a wide spectrum of SVD feature burden, stroke lesion size, and index stroke lesion type (i.e. cortical vs lacunar). The sample clinical characteristics have been published previously [16,17]; those relevant to this work are condensed in Table 1. The baseline hypertension (y/n) defined as a previous history of hypertension, or hypertension diagnosed at presented of stroke, age, and percentages of WMH in intracranial volume were retrieved  (1). Initially, we process each case by sampling all regions of interest using an anatomically-relevant template (2), describing local signal variations in regions of interest using descriptors (3). Subsequently, we study the dynamic descriptors using functional principal component analysis (4) and statistical analysis (5). ROI: region of interest. FPCA: functional principal component analysis. T 1 ,…,T TP : each time point of the DCE-MRIs. Subject 1 ,…,Subject P : each patient in the cohort. ROI 1 ,…,ROI R : each region of interest. from the study database. Additionally, we considered a visual clinical rating recorded at inclusion, total SVD [18] score, to account for four neuroimaging features of the SVD (lacunes, microbleeds, perivascular spaces, and WMH). Imaging was carried out on a 1.5 T MRI scanner (Signa HDxt, General Electric, Milwaukee, WI) using an eight-channel phased-array head coil. Both diagnostic and dynamic MR imaging acquisition parameters have been detailed in [19]. Diagnostic MRI at stroke presentation consisted of axial T2-weighted imaging (TR/TE = 6000/90 s, 24 × 24 cm field of view, 384 × 384 Propeller acquisition, 1.5 averages, 28 × 5 mm slices, 1 mm slice gap), axial fluid-attenuated inversion recovery imaging (TR/TE/TI = 9000/153/2200 ms, 24 × 24 cm field of view, 384 × 224 acquisition matrix, 28 × 5 mm slices, 1 mm slice gap), gradient echo imaging (TR/TE = 800/15 ms, 20°flip angle, 24 × 18 cm field of view, 384 × 168 acquisition matrix, 2 averages, 28 × 5 mm slices, 1 mm slice gap) and sagittal 3D T1weighted imaging (inversion recovery-prepared spoiled gradient echo TR/TE/TI = 7.3/2.9/500 ms, 8°flip angle, 330 × 214.5 cm field of view, 256 × 146 acquisition matrix, 100 × 1.8 mm slices). DCE-MRI acquisitions were obtained at approximately one month after first stroke presentation and consisted of a 3D T1w spoiled gradient echo sequence with TR = 8.24 ms, TE = 3.1 ms, 24 × 24 cm FOV, reconstruction matrix 256 × 192 and 42 × 4 mm slices. Following a precontrast acquisition, an intravenous bolus injection of 0.1 mmol/kg of gadoterate meglumine (Gd-DOTA, Dotarem, Guerbet, France) was administered with the start of 20 further acquisitions with 12°flip angle and a temporal resolution of 73 s, leading to a DCE-MRI duration of about 24 min (≈21 time points).

Image processing and region-of-interest sampling
We performed all image analysis blindly to clinical and permeability data. We aligned all time points of the DCE-MRI acquisition to the 12°p re-contrast image to correct for bulk patient movement using FSL-FLIRT [20,21]. For determining the WMH volume percentage in intracranial volume, we applied a segmentation method that has been evaluated previously against manual annotations in images acquired with similar scanning protocols [16]. On 150 subjects, the average difference on ICV was 2.7% (95%CI = 7%). On 20 individuals, the Jaccard similarity coefficient for WMH was 0.61 (95%CI = ± 0.37). In a test-retest analysis on 14 cases comprising volunteers and patients with mild non-disabling stroke, the coefficient of variation for repeated measurements of the segmentation technique was 0.21 [22]. Furthermore, trained analysts double-checked and manually edited these segmentation masks under the supervision of an experienced neuroradiologist.
We sampled five brain regions, comprising blood vessels [BL], CSF, deep GM [GMD], cortical GM [GMC], and WM, using circular nonoverlapping samples that covered approximately 12 mm 2 in-plane and were distributed throughout four slices [23][24][25], as described in Table 2 and exemplified in Fig. 2. The figure illustrates four slices and sampling points for one of the patients in our cohort, but it does not indicate sampling points are fixed nor predetermined for every single patient. In fact, these spots vary to avoid areas with partial volume effects, evident contrast-enhanced related truncation artefacts, WMH, enlarged perivascular spaces, mineral depositions, lacunes and ischaemic or haemorrhagic lesions to avoid biasing our analysis. The samples were initially placed by a trained analyst using Analyze 11.0 (AnalyzeDirect Inc., Mayo Clinic), edited by another one, and, finally, agreed between both observers. We opted for sampling brain regions to reduce the influence of the spatial densities, partial volume effects, and avoid obvious correlations/associations as a result of descriptors encoding volumetric information [26].

Power spectrum
We used the radial power spectrum (RPS) to account for spatial variations in the strength of the signal frequency components. First, we sampled the signal using the anatomical-relevant template. Second, we used the 3D discrete Fourier transform to obtain a representation of each volume in the frequency domain. Let I ∈ℝ N×N×N be a brain MR volume, the corresponding discrete Fourier transform, Third, we computed the magnitude spectra and averaged all the frequencies over concentric rings of width 1 using the following formula represent the corresponding spherical coordinates. For each time point and region of interest, signals were described using 256 frequencies.

Grey-level co-occurrence matrix based descriptors
We measured local signal variations using metrics of homogeneity and variability extracted from grey-level co-occurrence matrices (GLCM). These matrices summarise the co-appearance of intensity values in an image, i.e. they quantify the frequency at which two intensity values occur in the same neighbourhood. Mathematically speaking, a co-occurrence matrix is calculated as follows: where x ( ) denotes the set of voxels in the neighbourhood of x. The matrix C reveals information of the region of interest. For instance, the higher the values in the diagonal, the more homogeneous the region under examination.
Haralick et al. [27] proposed various measures of homogeneity and heterogeneity of region of interest based on the normalised GLCM values. The process for computing them was four-fold. First, we quantised each region of interest in 2 4 grey levels. Second, we computed the GLCM using an eight-connected neighbouring structure. Third, we normalised the GLCM by dividing each cell by the total number of voxel pairs. Fourth, we computed energy, contrast, correlation, variance, inverse difference moment, sum average, sum variance, sum entropy, and  (6), and WM (28) High Must be the first or second slice above the ventricles Superior sagittal sinus (1), GMC (6), and WM (28) entropy [27], as condensed in Table 3. For each time point and each region of interest, the number of GLCM features per patient was 9.

Local binary patterns descriptors
Local binary patterns (LBP) are another visual descriptor that quantify local signal variations. The original approximation characterises the way voxels relate to their neighbourhood using binary codes and summarises them for the entire region of interest using a histogram. The process is three-fold. First, the binary code for each voxel is computed by analysing the way its intensity relates to the ones of its neighbours: if the value is higher than the one of its neighbour, we assign a zero; and a one otherwise. For instance, if the voxel value is equal to 4 and the neighbouring intensities Second, the frequency of each one of these codes is totalised. Third, the resulting histogram is divided by its sum to express probability. The normalised histogram is used as a texture descriptor.
In this work, we considered two variants of the original approximation called uniform local binary patterns (ULBP) [28] and local configuration patterns (LCP) [29]. The former maps the original set of codes (2 8 = 256 different binary patterns) to a subset of 59 codes to reduce the cardinality of the histogram and provide the descriptor with a simple rotation invariance. The latter combines the LBP with rotation invariance and another descriptor which quantifies local linear dependencies between a voxel and its neighbours. We calculated these descriptors for each region of interest for each patient using a radius equal to 1. For each time point and each region of interest, the number of ULBP and LCP features per patient were 59 and 81, in that order.

Functional principal component analysis
Comparing the resulting dynamic descriptors could be an intricate task due to the dimensionality of the problem: Conventional feature reduction approaches, such as traditional principal component analysis (PCA) or auto-encoders, are not suitable for this problem since the data temporality would be neglected in principle. Therefore, we resorted to applying the functional principal component analysis (FPCA) method proposed by Happ and Greven [30] which allows us to model each descriptor as a function in time, reduce the time dimension respecting its nature, and obtain a single score per principal mode of variation. In a nutshell, the idea is to reduce the cardinality in one dimension (time) and then on another one (space). Each one of the elements of each descriptor can be seen as a function in time. In such a way, we could find the eigenvalues and eigenfunctions that better describe them. Let D be the number of elements under study, P = 42 the number of patients, and R = {R (1) ,…,R (D) } the set of elements, each of them described by the corresponding measurements, r j   Table 3 Considered texture descriptors based on the grey-level co-occurrence matrix. C a+b (k) represents the grey level sum distribution, expressed as C a+b (k) = ∑ a ∑ b {C(a,b) if a + b = k,0 otherwise}. SA in the sum variance formula denotes the sum average operation.

GLCM metric Brief description Formula
Contrast Overall contrast. The higher the co-occurrence of low and high intensity values, the higher the contrast.
Linear dependency between neighbouring voxels. Values close to 1 reflect high correlation.

Energy
Second angular moment measuring uniformity. Energy increases with increased homogeneity.
Dispersion of values around the mean. Variance increases with increased heterogeneity.
Measurement of randomness. Higher entropy values indicate higher heterogeneity.
Inverse difference moment Overall homogeneity. Higher values indicate higher homogeneity. ∑ Sum variance Dispersion of values around the sum average. Higher values of sum variance relate to higher heterogeneity.
. We set M (j) to five as resulting eigenvectors accounting for the 99% of the univariate variation. Third, all of these scores ξ ik j ( ) were arranged in a matrix form, ( ) . Fourth, scores were calculated using eigenanalysis on the covariance matrix of Ξ. We resorted to analysing the first mode of variation. The output was a single score per subject and per region of interest.

Validation against clinical parameters
We applied statistical tests to determine whether patients with similar health status exhibit similar principal component (PC) scores. We used the Kruskal-Wallis test for testing differences in PC scores between patients grouped by overall burden of SVD features. We considered multiple linear regression to establish whether age, WMH volume, biological sex, and stroke lesion type were associated with the PC scores (i.e. PC score = β 0 + β Age ⋅Age + β WMH ⋅WMH vol + β Sex ⋅Sex + β Lac ⋅Lac, where each β ∈ℝ is a standardised regression coefficient).

Comparison against other techniques
We compared our proposal against two surrogate measures of small vessel dysfunction: area under the enhancement curve (AUEC) [31] and signal enhancement slope [32,33]. The idea behind both approaches is as follows. The former computes integral of the enhancement curve over time since higher accumulation of contrast agent leads to higher enhancement (thus, higher AUEC). The latter calculates the slope of the enhancement curve (assuming linearity after the bolus arrival peak) relates to the degree of contrast agent leakage. In both cases, we computed the enhancement curve on the sampling points described in Section 2.2, measured AUEC and slope, and used the statistical tests in Section 2.5 to establish the relationship between these measurements and clinical variables.

Experiments and results
We applied our framework to the 42 cases by segmenting each of the 42 DCE-MRI scans, sampling the signal in each region of interest, describing signal in each time point and region of interest, and applying FPCA on the resulting dynamic descriptions to analyse directions of variations within the data. We explored whether these variations were associated with clinical variables.
We grouped the PC scores by the total SVD score and applied the Kruskal-Wallis test to determine whether there were significant differences between groups. The results are displayed in  = 4, respectively). The other four descriptors did not vary significantly with the overall burden of SVD in any of the tissues in any of the four slices: signal enhancement, GLCM metrics, area under the enhancement curve, and enhancement curve slope. We carried out multiple linear regression to investigate whether biological sex, age, WMH volume, and stroke lesion type were associated with the observed PC score. The results are condensed in Table 5. Overall, we observed that (i) these four covariates influenced the features in blood vessels, CSF, deep GM, and WM, but not the ones on the cortical GM region, and (ii) the analysis of the signal enhancement slope did not result in significant associations.
Regarding blood vessels (carotid arteries, basilar artery, sagittal sinus), the regression results implied strong relationships health status and observed intensity and spectral features (signal enhancement, RPS, and AUEC), explaining approximately 20% of the data variation and being significant predictors of these models (p-value < 0.05). Signal Table 4 Kruskal-Wallis values obtained by grouping PC scores by total SVD score. The results are expressed concerning χ 2 and Pr (p-value). The degrees of freedom were four. Enh, RPS, GLCM, ULBP, LCP, AUEC and slope stand for signal enhancement, radial power spectrum, grey co-occurrence matrix metrics, uniform local binary patterns, linear configuration model, area under enhancement curve, and enhancement curve slope, respectively. Values in bold are significant (p-value < 0.05).  (continued on next page) enhancement values were associated to the overall burden of WMH in low-middle (β = 0.31;p-value = 0.05) and high (β = 0.35;p-value < 0.05) slices and the biological sex of the patients in middle-high (β = −0.30;p-value < 0.05) and high (β = −0.29;p-value = 0.05) slices. These associations were similar for AUEC, except that the features measured in the low-middle slices were not significant (pvalue > 0.05).
In the low-middle slice, the overall burden of WMH contributed significantly to the observed spectral features in deep GM (β = 0.45;pvalue < 0.01) and WM (β = −0.37;p-value < 0.05), describing approximately 20% and 15% of their variability. However, no other descriptor seemed to encode any relevant information for these two tissues in this slice.

Discussion
In this paper, we propose a framework incorporating dynamic spectral and textural descriptors and functional principal component analysis to study dynamic brain MRI signals of brain pathology. In particular, we applied our processing pipeline to the study of SVD tissue changes using DCE-MRI acquisitions. We analysed the dynamic descriptors from blood vessels, CSF, grey and white matter brain regions of a population with features of SVD of differing extents to examine whether subjects with different biological sex, age, WMH volume, stroke lesion type, and overall load of SVD features exhibited distinctive patterns. To the best of our knowledge, this is the first time that these dynamic descriptors have been examined jointly for this purpose. Of note, our framework could be used for analysing other dynamic and non-dynamic brain MR acquisitions, with further testing.
Most of the features extracted from the CSF cavities near the choroid plexus (low-middle slice) were significantly associated with the biological sex or stroke lesion type, and the overall burden of WMH. We observed that the overall enhancement in the CSF cavities increased if patients were female or had a lacunar stroke, and with decreasing WMH burden. In accordance with previous findings [10], we observed a relationship between burden of features of SVD and leakage of Gadolinium-based contrast agent into CSF. However, our results suggest that the CSF enhancement might be inversely proportional to the burden of WMH. Bigger sample size and further testing are needed to determine the direction of the relationship. Our framework found features correlated with WMH burden and, as increased WMH load is known to be associated with blood-brain barrier leakage [9], our proposal shows promise for studying subtle small vessel dysfunction.
We observed that surrogate measures of small vessel disruption extracted from the signal enhancement curve, such as area under the enhancement curve, could capture relevant information linked with the health status of the patients, but some of the texture and spectral descriptors were more sensitive to variations in the deep GM or WM. In fact, we noticed that while both the slope and area under the enhancement curve measurements in WM did not reflect the health status significantly, but spectral and texture descriptors did. This might be a consequence of spectral and texture descriptors taking into account and encoding neighbouring relationships. The power spectrum reflected information associated with overall WMH burden in blood vessels, CSF, deep GM, and WM; being an encouraging outcome since it shows the analysis of the power spectra is more descriptive than current semiquantitative surrogates of dysfunction.
The current proposal exhibits two drawbacks: the region sampling strategy and the generalisability. First, region sampling prevents descriptors from encoding volumetric information, but sample selection is tedious and not resilient to motion. Due to the prolonged scanning process, brain DCE-MRI acquisitions are prone to head motion artefacts. In small vessel dysfunction assessments, imaging artefacts confound whether tenuous enhancements are a consequence of subtle blood-brain barrier abnormalities and, although the repercussions of motion have not been documented for this particular problem, we believe that they compromise both the interpretation and subsequent result interpretation. Notwithstanding that we realigned all time points to the precontrast scans to correct for bulk patient movement, this approximation does not compensate for possible information loss or k-space dephasing due to motion. Second, the imaging protocol influences the features that are captured by the different descriptors as none of them is scaleinvariant in principle and also the synthesis step as lower temporal resolution results in less information in the time domain. Even though the acquisition protocol was fixed in this study, their application to multi-centre studies might be restricted.
Future work should contemplate finding automatically anatomically-relevant regions proximal to arterial territories or in which signal to noise ratios are the highest. In such a way, we could enlarge our sample size to (i) draw stronger conclusions, (ii) establish clearer links between computer vision descriptors and underlying physiopathological processes, and (iii) determine whether these descriptors are useful in dysfunction classification problems. Additionally, the current proposal needs to be tested on diagnostic sequences and pre-vs post-contrast analyses to examine whether they capture dysfunction-related information.
Our findings add confidence to previous studies in which DCE-MRI signals from patients with different age, health status, and premorbid brain condition exhibited different tendencies [26,33]. Our proposed framework seems promising and feasible, but it needs further testing on a larger sample and on pre-vs post-contrast and cross-sectional studies.

Declaration of competing interest
None.