Spatial heterogeneity analysis of brain activation in fMRI.

In many brain diseases it can be qualitatively observed that spatial patterns in blood oxygenation level dependent (BOLD) activation maps appear more (diffusively) distributed than in healthy controls. However, measures that can quantitatively characterize this spatial distributiveness in individual subjects are lacking. In this study, we propose a number of spatial heterogeneity measures to characterize brain activation maps. The proposed methods focus on different aspects of heterogeneity, including the shape (compactness), complexity in the distribution of activated regions (fractal dimension and co-occurrence matrix), and gappiness between activated regions (lacunarity). To this end, functional MRI derived activation maps of a language and a motor task were obtained in language impaired children with (Rolandic) epilepsy and compared to age-matched healthy controls. Group analysis of the activation maps revealed no significant differences between patients and controls for both tasks. However, for the language task the activation maps in patients appeared more heterogeneous than in controls. Lacunarity was the best measure to discriminate activation patterns of patients from controls (sensitivity 74%, specificity 70%) and illustrates the increased irregularity of gaps between activated regions in patients. The combination of heterogeneity measures and a support vector machine approach yielded further increase in sensitivity and specificity to 78% and 80%, respectively. This illustrates that activation distributions in impaired brains can be complex and more heterogeneous than in normal brains and cannot be captured fully by a single quantity. In conclusion, heterogeneity analysis has potential to robustly characterize the increased distributiveness of brain activation in individual patients.

Introduction

Brain activation maps derived from functional MRI measurements are usually compared between diseased and control subject groups to draw statistical inferences on aberrant activation patterns related to the neurological disease conditions. Such inferences on aberrant brain activation patterns rely on the spatial overlap of activated regions among the members of a well-defined sub-population (Haxby et al., 2001; Manoach et al., 2000). However, even in the normal brain, cognitive functions depend on network activity of which the activation patterns are inherently distributed and are likely affected by natural variability of the brain’s functional organization, especially in more complex tasks. Any abnormalities in brain activation may not only express as aberrant activation levels but also as modifications in the distribution of activated regions. Inferences on abnormal brain activation are therefore not unambiguous as a priori we do not know whether differences in group averaged activation results are due to differences in activation level or reductions in overlap of activated regions. To provide further insight in aberrant activation maps of impaired brains, novel measures are required that characterize, and preferably quantify, the heterogeneity of activation patterns at the individual subject level. Such measures should ideally characterize the organization of the distribution of brain activation rather than the local degree of activation overlap.

Multiple measures for the characterization of spatial heterogeneity have previously been explored in other fields of medical imaging (Bright et al., 2009; Damona et al., 2008; Huettel et al., 2004). In the field of tumor imaging, spatial heterogeneity of contrast enhanced structural scans has for a long time been considered as a marker of malignancy (Alic et al., 2006; Alic et al., 2011; Jackson et al., 2007; Mohajen et al., 2011; O’Connor et al., 2011; Rose et al., 2009; Tixier et al., 2011).

Within the domain of fMRI, heterogeneity in time-series data has also been explored previously. Zang et al. (2004) presented a method for fMRI data analysis based on regional homogeneity, in which Kendall’s coefficient of concordance (Baumgartner et al., 1999) was used to compare the similarity of time-series in a voxel to the series from its neighbors. This method has also been applied in studies of Alzheimer's disease (Liu et al., 2008), neuromyelitis optica (Liang et al., 2011), aging (Wu et al., 2007), individual intelligence (Wanga et al., 2011), autism spectrum disorder (Shuklaa et al., 2010), depression (Yao et al., 2009), and attention-deficit hyperactivity disorder (ADHD). Leech and Leech (2011) quantified spatial heterogeneity in fMRI as a variation in the intensities of activated and neighboring voxels. These publications suggest that there could be increased activation map heterogeneity in certain patients. All the studies have investigated heterogeneity in activation patterns between subject groups, whereas this study aims to quantify spatial heterogeneity on the subject level.

Previous studies have tried to quantify the existence of heterogeneity in general. The current study shows how to characterize heterogeneity more specifically into different dimensions. We have evaluated well-known heterogeneity measures of functional MRI derived activation patterns in patients with epilepsy and language impairment and compared these to those obtained in age-matched healthy controls. The measures are applied to activation maps pertaining to language and motor tasks. The primary hypothesis is that the spatial heterogeneity of brain activation maps is stronger for patients than for controls. To further strengthen this hypothesis in relation to the disease characteristics, we secondarily investigated whether the increase in spatial heterogeneity is specific to the language-related, and not motor-related, activation pattern in (language-impaired) patients. The focus of the current study is on the application and evaluation of heterogeneity measures that can be applied to activation maps of individual subjects with a neurological disorder.

Methods

Overview

The selected heterogeneity measures focus on different aspects of spatial distributiveness, including the shape (compactness), complexity in the distribution of activated regions (fractal dimension and co-occurrence matrix), and gappiness between activated regions (lacunarity). The various heterogeneity measures are described below and are calculated for the activation maps of a language and a motor task. Furthermore, a support vector machine (SVM) classifier is used to analyze the discriminative power of the heterogeneity measures. The block diagram in Fig. 1 outlines the methodological steps used for the spatial heterogeneity analysis.

Fig. 1

Block diagram of the subsequent methodological steps in the heterogeneity analysis.

To exclude the influence of other factors, image noise and head motion analysis is applied to the activation maps. Head motion analysis is performed by computing displacement, entropy and smoothness of the three dimensional head movement time-series. The noise level is determined by the signal-to-noise-ratio (SNR) from the acquired echo-planar images.

Heterogeneity analysis

In this section we explain each method to determine the different heterogeneity measures in detail. In Appendix A the heterogeneity measures are illustrated using a conceptual set of images.

Overlap of activation

To obtain a measure for the degree of spatial activation overlap, we calculated the number of overlapping (i.e. commonly) activated voxels in the patient and control groups separately as a function of activation t-threshold. Note that overlap of activation represents a measure of heterogeneity at the group level rather than at the individual level.

Number of activated regions

The number of activated regions is compared between patients and controls as a function of the t-threshold, which controls the most significantly activated x-% (x varies from 3–10) of voxels.

Fractal dimension

A fractal dimension is a statistical index of complexity comparing how details in a pattern change with the spatial scale at which it is measured (Dougherty and Henebry, 2001). The concept of fractal dimension was first proposed by Mandelbrot (1983) to describe the geometry of objects in nature. The fractal dimension can be seen as a number that describe the texture (and thus heterogeneity) in the object. Objects with high fractal dimension are heterogeneous. There are various types of fractal dimension analysis techniques that could be used to describe the complexity (fractal pattern) of objects such as the frequently used box dimension, also called the Renyi dimension (Mandelbrot, 1983). The box dimension is a computationally simple method to quantify heterogeneity in binary images. It is computed by imposing regular grids of a range of spatial scales on the object studied and counting the number of grid elements (boxes) n that are occupied by the object. By plotting the number of occupied boxes as a function of the reciprocal scale on a log–log axis the box dimension represents the slope of the regression line.

Shape of activated regions

The shape of the activated regions is compared to a (single) perfectly spherical region. The measure that quantifies similarity of a region to a sphere is the isoperimetric quotient Q and is defined as

in which V and S are the volume and surface area of the activated regions, respectively. The volume is the total number of all activated voxels and the surface area is the number of activated voxels which are located at the surface of the activated regions. Voxels are classified to be on the surface of activated regions, if at least 1 out of their 26 neighboring voxels is not activated. The value of Q decreases when the region deviates from a sphere, because of the relative increase in surface area. Also the splitting of one large region into multiple smaller (distributed) regions, with preservation of total activation volume, will lead to a lower (averaged) Q. Note that the spatial distribution of multiple (non-overlapping) activated regions has no effect on Q.

Co-occurrence matrices

The co-occurrence matrix is defined as the distribution of co-occurring intensity values (t-values) at a given spatial offset (Amadasun and King, 1989; Thibault et al., 2009). This matrix is widely used for the characterization of texture in gray-scale images (Thibault et al., 2009). It is computed as

Here C is the co-occurrence matrix, n, m and o represent the sizes of the image and I(p,q,r) is the value at voxel position (p,q,r). The co-occurrence matrices used in this study are computed by taking mean of co-occurrence matrices computed at different offsets (Δx = 1, Δy = 0, Δz = 0), (Δx = 0, Δy = 1, Δz = 0), …, (Δx = D, Δy = D, Δz = D), where D is selected such that difference in entropy at D and D−1 is less than 0.0001. Two texture measures, i.e. homogeneity and entropy (Filho and Sobreira, 2008), were derived from these co-occurrence matrices. Homogeneity is the measure of image uniformity, which is inversely related to the heterogeneity, and is defined as

C(i,j) is the value of the co-occurrence matrix at the ith row and jth column. Entropy is a measure of randomness in an image and increases with heterogeneity and is defined as

Here C′(i,j) is the normalized value of the co-occurrence matrix at (i,j) by dividing each element of the co-occurrence matrix by the sum of all matrix elements C(i,j).

Lacunarity

Lacunarity is a measure of how sparse patterns fill space. It gives higher values where patterns have more, larger or more variable gaps (Dougherty and Henebry, 2001; Filho and Sobreira, 2008). Lacunarity was also originally developed by Mandelbrot (1983) to describe a property of fractals, and has since been extended to describe real data sets that may or may not resemble fractal distributions. Mandelbrot proved that different textures might have the same fractal values, but different lacunarity values. Therefore lacunarity can be interpreted as a complementary measure to the fractal dimension. Fractals depend on the geometric structure of an object whereas lacunarity is more sensitive to the sizes and distribution of gaps between the geometric structures. Geometric objects with low lacunarity are homogeneous because all gap sizes are the same, whereas high lacunarity objects are heterogeneous. It is important to note that objects that are homogeneous at a particular scale might be more heterogeneous at other scales, therefore lacunarity may depend on the spatial scale of measurement.

Lacunarity can be implemented in 3D using the sliding box approach. According to this algorithm a box of size r × r × r slides (overlapping) over a volume image. The mass M of an image is then defined by the number of activated voxels in a gliding box. n(M,r) represents the number of gliding-boxes with size r × r × r containing M activated voxels. The probability distribution Q(M,r) is obtained by dividing n(M,r) by the total number of boxes B:

Lacunarity at scale r is defined as the mean-square deviation (i.e. second moment) of the variation of mass distribution probability Q(M,r) divided by its square mean:

Conceptually, this resembles the (relative) variability of activated region sizes normalized to the (squared) average region size at a certain spatial scale r. We have used multiple cubic gliding box sizes ranging from 8 × 8 × 8 mm to 16 × 16 × 16 mm.

Intelligent parameter selection and classification

Heterogeneity is a complex feature that may not be captured by just one measure but needs multiple measures. To identify and combine the most relevant heterogeneity measures, a classification method was applied. Following Chen and Lin (2005), the F-score has been used to compute the importance of heterogeneity measures for distinguishing patients and controls. It is defined as

Here, xk, k = 1, …, T, is a training vector and T is total number of samples. P and C are the total number of patients and controls, respectively. are the mean of the ith measure of the total, patient, and control data sets, respectively; , is the ith measure of the kth patient, and , is the ith measure of the kth control. The numerator indicates the discrimination between the patients and controls (interclass variation), and the denominator indicates the discrimination within each of the two sets (intraclass variation). Higher F-scores indicate better discrimination.

Given a set of parameters for each subject, our objective is to classify each subject as a patient or control based on the values of the measures. To perform the classification based on the calculated heterogeneity measures, we have used the support vector machine (SVM) classification algorithm (Kumar, 2004; Zeng et al., 2007). SVM has been applied with a polynomial kernel (degree 2), using the selected measures with the F-score method. A significant advantage of the SVM algorithm is the search for a global and unique solution.

Statistical analysis

Differences in heterogeneity measures between patients and controls were assessed using a non-parametric Wilcoxon rank test (Conover, 1980). A non-parametric test has been used because distribution of the data is unknown. To determine the degree to which extent the different heterogeneity measures were independent of each other, the (Pearson) correlation matrix over all measures was calculated.

Furthermore, the sensitivity and specificity of the different heterogeneity measures were calculated to determine the potential of the different measures to distinguish the activation maps of patients from controls. This analysis was performed using 5-fold cross validation (Kohavi, 1995). For each fold, all the data were randomly assigned to five sets S1, S2, …, and S5 so that all sets were of equal size and each set contains an equal number of patients and controls. Then training was performed on S1, S2, S3, and S4 and testing on S5, followed by other combinations such that each set has been used once for testing. This approach has the advantage that the training and testing sets were both relatively large, and each data set is used for both training and validation on each fold.

Receiver–operator-characteristic (ROC) analysis of all heterogeneity measures was used to determine the discriminative power of each heterogeneity measure. For each heterogeneity measure, the area-under-curve (AUC) was calculated by varying the thresholds. The AUC was also obtained using 5-fold cross validation to provide the uncertainty.

The spatial correspondence between the activation maps of the two groups was characterized by the Jaccard index, which is the measure of similarity. This index was calculated as a ratio of the number of overlapping activated voxels to the total number of distinct voxels. Jaccard index will be lower for the group with less overlapping voxels, i.e. group with more spatial variability in activation.

Motion artifacts

To test whether patients give rise to stronger motion induced confounds in functional images, for instance spurious activation regions, the absolute and relative head motion time-curves were evaluated in terms of displacement, smoothness and entropy and correlated to the activation heterogeneity measures. The description and results of this analysis are given in Appendix B.

Noise analysis

SNR was determined on the original T2*-weighted images for each subject. The region inside the brain is considered as the signal and outside the brain as the noise. SNR is computed by dividing mean value signal inside the brain by the standard deviation of the signal outside the brain.

Subjects

Twenty-three children were recruited from the database of a specialized tertiary referral center for epilepsy upon having a clinical diagnosis of Rolandic epilepsy. Rolandic epilepsy is the most common benign childhood epilepsy, has a genetic basis and is characterized by centrotemporal spikes on the electroencephalogram (EEG) (Panayiotopoulos et al., 2008). Seizures occur mostly during the night and involve hemifacial spasms and speech arrest. In recent years evidence has accumulated that Rolandic epilepsy is associated with co-morbidities such as language impairment, despite its mild seizure semiology.

The average age at testing was 11.4 years (range 8–14 years) and we included 7 girls and 16 boys. A control group of 21 age and gender matched healthy children was also included (10 girls, 11 boys, mean age 10.3 years, range 8–14 years).

A previous study on these groups revealed a bilateral activation pattern with no clear (overlapping) differences in activation maps between the children with Rolandic epilepsy and healthy peers (Besseling et al., 2013).

Language assessment

The Clinical Evaluation of Language Fundamentals (CELF-4) test for children, Dutch edition (Paslawski, 2005; Semel et al., 2010), was used to assess language performance (Overvliet et al., 2013). The central outcome measure of the CELF test is the core language score, which can serve as a screening measure for language impairment. The patients' core language score was 92 ± 18(mean ± SD), which is significantly lower than the norm value of 100 ± 15(p = 0.047).

MRI

For fMRI a blood oxygenation level dependent (BOLD) sequence (T2*-weighted) was used at 3 Tesla field strength with the following settings: single-shot echo planar imaging (EPI), echo time/repetition time (TE/TR) 35/2000 ms, pixel size 2 × 2 mm2 and 4-mm thick axial slices. Each scan consisted of 195 whole-cerebrum dynamic acquisitions. The acquisition time was 6.5 min.

Language paradigm

A standard block design was used consisting of 6 task condition blocks interleaved with baseline condition blocks. Each block lasted 30 s and each paradigm started and ended with a baseline block.

The task comprised a word generation paradigm. Subjects had to covertly generate as many words as possible starting with a visually presented letter (U–N–K–A–E–P). The paradigm consisted of 6 word generation blocks (1 letter per 30-s block) alternated with baseline rest blocks (30-s), in which an asterisk (*) was presented for eye fixation. Previous studies have demonstrated activation in the anterior cingulated and inferior and middle prefrontal cortex in adults (Backes et al., 2005; Deblaere et al., 2002; Price, 2010).

Motor paradigm

In the motor task, left handed finger tapping was alternated with right handed finger tapping. An arrow was visually presented to indicate which hand to use. Both tapping conditions were interleaved with baseline rest blocks (each block lasts 25 s) and activation maps were derived by contrasting the tapping conditions against baseline.

Functional image processing

All fMRI time series were processed to generate activation t-maps using the Statistical Parametric Mapping (SPM8) (Friston et al., 1994) software application (Welcome Department of Cognitive Neurology, London, UK), which consisted of several steps.

First, the images of each dynamic series were realigned to correct for head motion using rigid transformations (translation and rotation, 6 degrees of freedom); these motion parameters were compared between groups. Next, the T1-weighted anatomical reference image was segmented into gray matter, white matter and cerebrospinal fluid components and transformed to the Montreal Neurological Institute (MNI) standardized stereotactic coordinate system. This normalization was subsequently applied to the realigned dynamic images. To correct for spatial misregistrations and to strengthen the assumption of normal distribution of data, the dynamic images were spatially smoothed using a Gaussian kernel with a full width at half maximum (FWHM) of 6 mm.

The task design was convolved with a standard hemodynamic response function to model the fMRI time series. Activation maps were statistically evaluated using t-contrasts.

Second-level analysis was performed to test for differences in overlapping task responses between the two subject groups.

Results

The values of the heterogeneity measures are listed in Table 1 for both tasks. For the language task, the number of activated regions, fractal dimension, entropy and homogeneity are slightly higher for patients than controls, but the differences were not significant. The isoperimetric quotient provided significantly lower values for patients compared to controls. Lacunarity showed significantly higher values for patients than controls. For the motor task none of the heterogeneity measures showed a statistically significant difference between patients and controls. Fig. 2(a)–(f) shows the heterogeneity measures of patients and controls as a function of percentage of most significantly activated voxels. All measures showed that the difference between patients and controls is conserved and comparable over activation thresholds. Fig. 3(a)–(f) shows the heterogeneity measures of patients and controls with t-values used as threshold. For increasing a t-threshold less voxels are activated. In Figs. 2 and 3 shape and homogeneity values decrease with the decreasing numbers of activated voxels, whereas lacunarity, entropy and fractals increase.

Table 1

Heterogeneity measures for pediatric patients with epilepsy and healthy controls computed with voxels top 5% of t-values, sensitivity and specificity of distinguishing patients and controls, the area-under-curve of the receiver–operator characteristic and the F-score.

Method		No. of regions	Fractal dimension	Isoperimetric quotient	Entropy	Homogeneity	Lacunarity	SVM
Language task	Patients	44.8 ± 2.9	2.45 ± 0.001	0.108 ± 0.0001	0.280 ± 0.003	0.973 ± 0.001	3.22 ± 0.22
		Controls	38.5 ± 3.6	2.45 ± 0.001	0.111 ± 0.0001	0.287 ± 0.002	0.974 ± 0.001	2.60 ± 0.13
		p-value	0.14	0.79	0.08	0.41	0.50	0.03
	Sn	78.2	86.9	86.9	21.7	34.7	73.9	78.3
	Sp	45.0	25.0	45.0	80.0	80.0	70.0	80.0
		AUC	0.59 ± 0.01	0.45 ± 0.01	0.63 ± 0.02	0.42 ± 0.01	0.56 ± 0.01	0.69 ± 0.02
		F-score	0.03	0.01	0.10	0.04	0.01	0.12
Motor task	Patients	50.8 ± 2.1	2.90 ± 0.001	0.0681 ± 0.001	0.63 ± 0.006	0.94 ± 0.001	10.7 ± 0.5
		Controls	54.5 ± 3.5	2.90 ± 0.001	0.0695 ± 0.001	0.63 ± 0.006	0.94 ± 0.001	11.3 ± 0.6
		p-value	0.42	0.55	0.91	0.60	0.47	0.56

Fig. 2

Heterogeneity measures in pediatric patients with epilepsy and healthy controls as a function of the percentage of the most strongly activated voxels (3–10%). For the fractal dimension (b), lacunarity (f) and entropy (d), the heterogeneity values increase with higher voxel percentages as the complexity, gap size and the distributiveness in the image increase. These measures are higher in patients than in controls at all the voxel percentages. The isoperimetric quotient (c) and homogeneity (e) decrease for higher voxel percentages as the activated regions become more irregular and distributed. It can also be observed that over a wide range of thresholds, activated regions in controls always appear closer to the ideal spherical shape and more homogenous than in patients.

Fig. 3

Heterogeneity measures in pediatric patients with epilepsy and healthy controls with t-values used as threshold on number of voxels activated. Similar to Fig. 2, for the fractal dimension (b), lacunarity (f) and entropy (d), the heterogeneity value increases with the increase in number of activated voxels and the isoperimetric quotient (c) and homogeneity (e) decrease.

The correlations between the different heterogeneity measures are listed in Table 2. All measures, other than the number of activation regions, show moderate to high correlation (0.7–0.9). Strong negative correlation was found between entropy and homogeneity (both based on the co-occurrence matrix).

Table 2

Table showing the correlation between different measures used in the study for language tasks. It can be noted that other than number of activated regions, all the measures show moderate to high correlation.

AUC values for the ROC of each measure are also given in Table 1. The highest AUC was obtained for the lacunarity analysis. The lacunarity measure also performed best with a sensitivity of 74% and a specificity of 70%. Using the SVM classification, the sensitivity increased to 78% and the specificity to 80%. The best two measures were lacunarity (F = 0.12) and isoperimetric quotient (F = 0.10). The F-scores of the heterogeneity measures are also given in Table 1. The highest 1% t-values are concatenated to form a feature vector, which was used for the SVM classification. The sensitivity and specificity achieved using t-values as features were 45% and 70%, respectively, which further indicates the importance of using heterogeneity measures.

Noise, motion and overlap analysis

The SNR of the patient group (36.4 ± 0.9) was not significantly different from the control group (34.9 ± 1.0).

Fig. 4(a)-(b) shows the number of clusters and mean cluster size as a function of the percentage of the most strongly activated voxels in patients and controls. The two figures illustrate that the number of clusters is higher in patients than in controls for a range of t-values, whereas mean cluster size is lower. To illustrate the distribution of activation levels, in Fig. 4(c) the histogram of activation maps for a range of t-values is shown. As these distributions do not differ, the total amount of activation is comparable between patients and controls.

Fig. 4

(a) Number of clusters; (b) mean cluster size for language task as a function of the percentage of the most strongly activated voxels in patients and controls and (c) histogram of t-values. (a,b) Illustrate that number of clusters are more in patients than in controls for a range of t-values whereas mean cluster is less in patients than in control. This effect shows that patients' activation pattern is more heterogeneous (distributed) than controls, however the effect is not statistically significant using only the number of clusters. (c) Shows that statistically the difference in number of activated voxels between patients and controls at different intervals of t-values is insignificant.

Fig. 5(a)–(d) shows the average activation maps of the patient and control groups for the language and motor tasks, which were highly similar. The group comparison yielded no significant differences. Fig. 6 shows a group comparison of the Jaccard index between patients and controls. The Jaccard index in patients appears to be lower than in controls indicating that the number of overlapping voxels in patients is lower than in controls and the distribution of activated voxels in patients is more heterogeneous.

Fig. 5

Activation (overlap) maps of the patient (a,c) and control (b,d) groups for the language (a,b) and motor (c,d) tasks. Depicted are the significantly activated voxels in color (p < 0.05). Similar regions are activated for the language and motor tasks in patients and controls. Second level analysis did not provide any statistically significant differences between the groups for both tasks.

Fig. 6

Jaccard index of overlapping voxels in patients and controls as a function of activation threshold.

Discussion

Current findings

Different well-known spatial heterogeneity measures were explored to capture the complexity, shape and spatial distribution of activated brain regions on an individual subject basis. The rationale of the current work was to determine spatial heterogeneity measures that can quantify the heterogeneity of brain activation. We evaluated the proposed measures in patients with epilepsy and language impairment and healthy controls in response to a language and motor task. All proposed measures of spatial heterogeneity revealed that patients exhibit a more heterogeneous activation pattern than controls in response to the language task. Shape- (i.e. isoperimetric quotient) and lacunarity-based measures performed better than other methods in distinguishing activation patterns of patients from controls. Similar heterogeneity analyses were performed for the motor task data. However, for the motor task heterogeneity measures did not reveal any differences. This indicates that differences in heterogeneity measures are most likely due to the responses of the impaired language function. The increased heterogeneity was not due to differences in image noise or head motion confounds.

The Jaccard index in patients was lower than in controls, and the spatial distribution of activated voxels in patients is more heterogeneous than in controls.

The current study provides neuronal correlates, in addition to neuropsychometric measures, for the functional deficits of the epilepsy disorder. A goal of such a neuroimaging study is that it adds to the understanding how the brain tissue responds differently in children with epilepsy with language comorbidity in comparison to healthy controls. As a clinical perspective, this may also serve in future clinical trials as a neuronal biomarker to support any treatment and development of the language co-morbidity.

Group averaged activation patterns

The degree of activation map overlap in the patient group was lower than in the control group. It was also seen that the number of activated regions in patients was higher than in controls, despite the observation that standard (overlap based) group comparison revealed no differences. These observations suggest that spatial heterogeneity in patients is stronger than in controls. However, the quantitative results of the overlap and the number of activated regions appeared not significant. Moreover, this overlap measure provides a group-based measure, which may describe strong variations between subjects of a group but not features of spatial heterogeneity in individual subjects. To overcome this limitation, we explored measures of spatial heterogeneity that can be applied to individual subjects.

Nature of spatial heterogeneity in brain activation

The two most important features in spatial heterogeneity differences between individual patients and controls were variations in shape (isoperimetric quotient) and sparsity (lacunarity) of activated regions. As the isoperimetric quotient was lower in patients, the distribution of brain activation regions in patients reflects activated regions with a surface area that deviates more strongly from a compact spherical region than for healthy controls. As the volume of brain activation was regulated by selecting consistently the highest t-value voxels, the shapes of the distributed activated regions are less compact (more irregularly shaped and smaller regions) in patients than in controls. The observation that the lacunarity measure appears elevated in patients demonstrates that activated regions are more separated and more irregular than in healthy controls. The number of activated regions as well as co-occurrence matrix derived measures (i.e. homogeneity and entropy) of spatial heterogeneity were less sensitive to differentiate the activation patterns of patients from controls. Co-occurrence matrices measure the texture content, in particular the repetition of a pattern in an image. The activation patterns in both patients and controls are more close to a stochastic texture pattern; therefore co-occurrence matrix measures do not provide a strong distinction between patients and controls.

Most remarkable is the observation that activation maps were more heterogeneous for the patient group than for the controls for the language task. Apparently, the increased heterogeneity of activation distribution was related to the impaired language function as for the motor task no differences in heterogeneity were found. However, more research with more tasks is required to draw final conclusions on this matter.

Dependence of the different heterogeneity measures

The combination of measures provides a better discrimination between activation patterns of patients and controls than using a single heterogeneity measure. The combination of measures indicates that the spatial organization of the activated brain regions is a complex feature that cannot be captured completely by one single heterogeneity measure. This is also observed in the correlation matrix for all the measures. All measures, other than the number of activated regions, were moderately correlated because they all capture to some extent the distributiveness in the activation map. However, the correlation is not perfect, which also points out that one single measure is not sufficient to capture the heterogeneity in patients and controls. This indicates that spatial heterogeneity in fMRI activation maps is due to the variation in all the three measures, i.e. shape, gappiness and complexity of the activation regions.

Interpretation of increased spatial heterogeneity

The stronger spatial heterogeneity of activated regions in patients compared to controls is observed in shape and spatial distribution of activated regions. This can be explained in the light of brain function. The brain is organized in functional networks of specialized cortical and subcortical regions. One of the interpretations of increased spatial heterogeneity in patients is that auxiliary brain regions may co-activate to compensate for disturbances in specific functional networks (Dupont et al., 2001; Eliassen et al., 2008; Vlooswijk et al., 2010). Such a compensation mechanism creates more irregular and smaller regions in the activation maps and also increases spatial dispersion. Another explanation for the increase in heterogeneity is that some areas in the impaired brains could activate even when that part of the brain is commonly activated for task performance (Eliassen et al., 2008). This may affect (suppress) the usual part of the brain to activate, and may increase the spatial heterogeneity.

The selection of a predefined percentage of the strongest activated voxels might include more false-positives that are randomly distributed throughout the brain. However, analyzing the heterogeneity measures by thresholding on activation t-value rather than percentage of activated voxels provided similar results. Moreover, the distribution of t-values, the image noise and head motion parameters were comparable for the patient and control groups.

Conclusion

Novel methods to quantify spatial heterogeneity in brain activation maps were explored. These methods focus on different aspects of heterogeneity of the distribution of activated regions, such as shape, complexity, geometric structure and gaps, which can be employed at the individual subject level. An increase in spatial heterogeneity was observed in pediatric patients with epilepsy relative to healthy controls, whereas no significant activation map differences were observed when employing conventional (overlap based) group comparison. Spatial heterogeneity in patients was best quantified in terms of shape and lacunarity measures, describing the distribution of the increased surface areas and the irregular gaps of the distributed activation pattern, respectively. We propose that spatial heterogeneity is a valuable functional (neuroimaging) biomarker that could be added to regular brain activation map analysis to characterize conditions of cerebral disease and should further be explored for other neurological disorders as well. Such neuronal biomarkers may also help in future trials that aim to diagnose and monitor patients that do or not respond to therapeutic assessments that aim to improve the language and other cognitive performances. We also expect that the applicability of the heterogeneity measures can be easily extended to other neuroimaging methods, for instance perfusion and diffusion imaging.

Table B.1

Displacement, smoothness and entropy parameters of the patient and control groups for the language task. None of the comparisons appeared statistically significant.

	Displacement		Smoothness		Entropy
	Absolute (mm)	Relative (mm)	Absolute (s−1)	Relative (s−1)	Absolute	Relative
Patient	0.57 ± 0.34	0.13 ± 0.16	3.2 ± 2.2	6.4 ± 3.6	7.32 ± 0.12	6.99 ± 0.45
Control	0.65 ± 0.38	0.13 ± 0.17	2.9 ± 1.3	6.9 ± 3.8	7.42 ± 0.10	7.12 ± 0.27

Table B.2

Correlation between heterogeneity measures and displacement, smoothness and entropy. All heterogeneity measures poorly correlate with displacement, smoothness and entropy. These measures are computed on the language task.

	Displacement		Smoothness		Entropy
	Absolute	Relative	Absolute	Relative	Absolute	Relative
Number of regions	0.27	0.31	0.32	0.21	–0.20	–0.31
Fractal dimension	0.12	0.16	0.11	–0.09	–0.17	–0.18
Shape	–0.04	–0.18	–0.02	0.19	0.19	0.14
Entropy	0.10	0.15	0.09	–0.02	–0.24	–0.23
Homogeneity	–0.08	–0.10	–0.07	0.01	0.22	0.22
Lacunarity	0.01	0.3	0.02	–0.26	–0.17	–0.04