Utilizing anatomical information for signal detection in functional magnetic resonance imaging

2.1 Linear models for voxel-wise multiple tests

Let Y_ixt denote the observed data from a functional MRI experiment at voxel $x$ and time $t$ for the $i$ -th subject. Here, we adopt the common view (cf. Section 5.4 in Ref. 28) of a mass-univariate linear model

The SPM then forms a random $t$ -field (cf. Ref. 53) with an inherent multiple comparison problem due to the large number of local hypotheses. One common strategy is to define local p-values at each voxel $x$ and for each subject $i$ based on the local values of the random $t$ -field and to control the family-wise error rate (FWER) using accordingly adjusted thresholds; cf. Ref. 54. However, this is known to be a very conservative approach with respect to the detectability of significant brain signals in the outlined framework. In contrast, approaches related to the control of the false discovery rate (FDR) can handle the multiple comparison problem, e.g., by the procedure proposed in. Ref. 6.

2.2 Parcellation of the human brain

Neuroanatomic research has found that the human brain can be parcellated into different sub-regions based on structural similarities. One of the earliest atlases is the Brodmann atlas (see⁷) which is based on the cytoarchitectural organization of the brain. The Brodmann Areas have been schematically transferred to a template brain, the so-called Talairach-Atlas (see Ref. 46) which is commonly used in fMRI studies to report the location of significant grand average activation, as used in Refs. 42, 43. For the analysis of these data in the current paper, we chose the Harvard-Oxford brain atlas (cf. Refs. 10, 31) that provides a parcellation based on gross anatomical landmarks and delivers an Aparc label $j$ for each voxel of each individual brain space. However, any other brain parcellation to define regional labels could be used with our methodology.

2.3 Statistical inference

As outlined in the introduction, we re-used fMRI data from a program code comprehension task first analyzed in Ref. 43 and performed a new analysis comprising the four steps outlined below. The experiment used two different levels of software program code comprehension stimuli, henceforth denote as bottom-up and top-down comprehension, to infer on the related cognitive processes. In our strategy, we combined the methods from Ref. 43 (steps 1 and 2) and Ref. 40 (step 3). Furthermore, we implemented our new methodological contribution of combining the evidence for activation of a given brain region across the subjects (step 4). For the first two steps, we conducted, for each subject, a random-effects linear model analysis as described above for deriving voxel-wise p-values.

Step 1: Program comprehension versus rest

In the first step, we contrasted (for each participant separately) the comprehension of program code (cf. Ref. 43) to the rest condition. This identifies brain areas with a positive deflection of the BOLD response. Furthermore, in order to account for the multiple comparison problem, we performed the Benjamini-Hochberg test (see Ref. 6) for FDR control. Only those voxels which have been declared significant by this procedure were considered in step 2. This methodology is justified by the fact that the FDR is an established screening criterion for high-dimensional multiple test problems.

Step 2: Bottom-up comprehension versus control condition

In this step, we contrasted (again, for each participant separately) one type of program comprehension: bottom-up comprehension. Bottom-up comprehension is induced when program code provides no semantic cues and programmers need to comprehend each line separately and then integrate the information in a slow, tedious process. For the significant voxels from the first step, we applied in a second step the same multiple test to the contrast of bottom-up comprehension against the control condition (syntax task) on the restricted set of voxels. As a result of this step, we get for each participant $i$ and for each considered voxel $k$ a p-value ${\tilde{p}}_{\mathit{ik}}$ .

Step 3: Regional p-values for every participant $\mathit{i}$

This step builds upon the methodology from Ref. 40, and it delivers for each participant $i$ and for each anatomical (regional) label $j$ a (confirmatory) significance evaluation with respect to the contrast specified in step 2. Hence, the evidence from all voxels of participant i in the brain region labeled by $j$ is combined in this step of data analysis.

To this end, let $\kappa$ be a tuning parameter with values in the interval [0, 1] (i.e. in per cent) and let $m_j$ be the number of voxels contained in the brain region labeled by region label $j$ . To keep the notation feasible, we implicitly assume here that for each participant $i$ the same number $m_j$ of voxels belong to the brain region labeled by $j$ . We consider the null hypothesis $H_{\mathit{ij}}$ of no relevant differential activation of the region labeled by $j$ for participant $i$ during the two tasks mentioned in step 2, together with its two-sided alternative hypothesis $K_{\mathit{ij}}$ . We call $H_{\mathit{ij}}$ the “regional null hypothesis” for the brain region labeled by $j$ for participant $i$ . We formalize $H_{\mathit{ij}}$ as a so-called partial conjunction hypothesis (see⁴⁰ and the references therein for a formal mathematical description), meaning that we consider the differential activation in region $j$ for participant $i$ relevant, if it contains at least $u_j≔\kappa ·m_j$ significant voxels. For testing $H_{\mathit{ij}}$ we calculate the “regional p-value” $p_{\mathit{ij}}^{\mathit{\text{REGION}}}$ , given by

In order to achieve family-wise error rate (FWER) control, we have to choose the tuning parameter $\kappa$ smaller than or equal to $1/J$ , where $J$ is the number of regional labels. Choosing $\kappa =1/J$ corresponds to the so-called Bonferroni multiplicity correction. The choice of $\kappa$ is discussed further in Appendix S1 of Ref. 40.

Step 4: Combined regional hypothesis tests by Fisher’s method

In this final step, we combine for each regional label $j$ the regional p-values calculated in step 3 over all participants $i=1,\dots ,n$ . In order to do this, we apply the so-called Fisher method to combine p-values. Namely, the Fisher test statistic $T_j$ for region $j$ is given by

Under independence of the data with respect to the participants, $T_j$ is asymptotically $X_{2n}^2$ -distributed (chi squared) with $2n$ degrees of freedom under the null. The latter independence assumption is justified, because the participants have been included in the study independently from each other.

Finally, we can reject the (over all participants $i$ combined) regional hypothesis $H_j$ (i.e., the respective partial conjunction hypothesis, but now with respect to the population, not with respect to a single participant) if and only if Fisher’s test statistic $T_j$ is larger than the (1-ακ)-quantile of the $X_{2n}^2$ -distribution with $2n$ degrees of freedom, where the tuning parameter $\kappa$ has been introduced in step 3. This parameter addresses the multiplicity of the test problem with respect to the $J$ regional labels which are simultaneously under consideration.

3.1 Simulation setting

We have simulated a dataset $\left\{Y_{\mathit{ixt}}\right\}$ with eleven subjects (referring to the index $i$ and corresponding to the group size in the real dataset below), a spatial grid of size $20\times 20\times 20=\mathrm{8,000}$ voxels (referring to the index $x$ ), and 195 time points (scans, referring to the index $t$ ). We have assumed that the 8,000 voxels are grouped into eight anatomical regions of size $10\times 10\times 10=\mathrm{1,000}$ each, and that these regions are correctly annotated for all eleven subjects. Two alternating stimuli in an ON-OFF block task design with a total of six ON blocks of a duration of 15 scans have been used for the temporal signal: In a circular area in one of the predefined anatomical regions the signal of one stimulus was twice as high as the signal by the other stimulus mimicking a signal contrast. In a corresponding area in a second anatomical region the signal was created with no difference between the stimuli. The datasets with first order autocorrelated Rician noise have been created using the R package neuRosim⁵¹; p-values ${\tilde{p}}_{\mathit{ik}}$ for Contrast 1 (ON of either stimulus versus rest) and Contrast 2 (one stimulus versus the other) where determined using the R package fmri⁴⁵. This simulation setup has been run 500 times.

3.2 Considered data analysis methods

For the statistical analysis of the simulated data, we have considered three different methods.

3.3 Results

We have assessed the type I and the type II error behavior of the three methods described in the previous subsection by calculating (for each of the three methods) the proportion of simulation runs in which any region without true activation has been declared significantly activated (type I error component) as well as the proportion of simulation runs in which the one region with true activation has not been declared significantly activated (type II error component). Table 2 summarizes the results of our computer simulations.

From Table 2 it becomes apparent that all three methods are capable of protecting reliably against type I errors under our data-generating model. However, the power of the proposed method appears to be considerably larger than that of the voxel-wise method. Namely, during our 500 simulation runs, our new method as well as the cluster-wise method always detected the truly activated region, while the voxel-wise method detected it only in 367 out of the 500 runs.

4.1 Previous findings

In the two previous studies mentioned before, Siegmund et al. used similar analysis processes. They used BrainVoyager™ QX 2.8.4[1]. The anatomical scans were transformed into the Talairach brain to account for differences in brain size; cf. Ref. 46. They preprocessed the functional data of both studies with a standard pipeline of: 3-D motion correction, slice-scan-time correction, and temporal filtering. In addition, they applied a spatial smoothing with a Gaussian filter (FWHM = 4 mm).

In the first study⁴², the random-effects GLM revealed five brain areas (BAs 6, 21, 40, 44, 47) with significant activation with the contrast Bottom-Up Comprehension versus Control condition, i.e. syntax task. In the second study, the same contrast revealed no significant areas anymore, likely due to the reduced statistical power of five instead of twelve bottom-up comprehension tasks per session. Thus, the authors ran a regions-of-interest analysis restricted to the identified activation clusters of the first study on the data of the second study. This resulted in a significantly stronger activation for Bottom-Up Comprehension versus Syntax in BAs 21, 40, and 44.

4.2 Data export and preparation for re-analysis

In Figure 1, we illustrate the overall process for the re-analysis of the data from Ref. 43. For our re-analysis, we exported the already pre-processed data from the second study. We did not use the data of the first study as our method does not need prior definitions of regions of interest. We used BrainVoyager to access the data from Siegmund et al. We exported the statistical values (i.e., t-scores, p-values) of the obtained brain activation on a voxel basis for each participant. The voxel resolution is the same BrainVoyager uses, i.e., a 1 mm interpolated resolution.

In addition, we used FreeSurfer to segment and parcellate the brain of each participant based on their anatomical scan; cf. Refs. 17, 18. We used the Destrieux’ cortical atlas to assign Aparc labels on an individual participant basis (see Ref. 11) as a suitable region definition. Next, we used Nipype (see Refs. 16, 21) to convert Freesurfer labels to a BrainVoyager-readable format.

Our last step annotated the exported functional data with the individual anatomical labels for each participant. We removed all functional voxels for coordinates that had no assigned Aparc label, which typically are voxels that are not considered gray matter.

4.3 Results

Figure 2 displays the six brain regions which have been declared as significant (at FWER level $\alpha =5\%$ ) associated with the task at hand by our described methodology. There is an overlap with the results obtained by Ref. 43 in which they utilized prior knowledge, but also some differences which we discuss in the following in their nuances. We visualize the confirmed network of brain activation from Siegmund et al. in Figure 3 and our identified network of significantly activated Aparc labels in Figure 4. Figure 3 and Figure 4 show overlapping results with regard to the Brodmann area 21 that covers the middle and inferior temporal gyrus (separated by the inferior temporal sulcus). However, we observed differences regarding smaller brain regions. We compare Siegmund’s replication efforts to our results in Section 5.1.

Figure 2. The six significant brain regions (at FWER level α = 5%).

In each row, each point corresponds to the Aparc p-value $p_{\mathit{ij}}^{\mathit{\text{APARC}}}$ of one study participant i, where the index j refers to the area indicated by the code at the beginning of the row.

Figure 3. Network of left-lateralized confirmed brain areas. In Ref. 43, BAs 21, 40, 44 were found activated during program comprehension.

Figure 4. Results of our analysis with significantly activated Aparc labels.

Activation are particularly in the middle and inferior temporal lobe.

Our method found three Aparc labels in the inferior and middle temporal gyrus in the left hemisphere. Siegmund et al. found their largest and most robust activation cluster in BA21 of the left hemisphere, which covers several gyri in the temporal lobe. These left temporal gyri are often associated with semantic processing of natural language, which is typically left-lateralized for right-handed participants. In the context of programming, the activation is believed to be responsible for extracting the meaning of individual identifiers and symbols during program comprehension. We found two further Aparc labels bilaterally in the inferior temporal gyrus and the anterior collateral sulci, which are both in the temporal lobe as well.

For each of these six regions indexed by $j$ , we display all subject-specific regional p-values $\{p_{\mathit{ij}}^{\mathit{\text{REGION}}}:1\le i\le n$ }, where $n$ is the number of study participants. For each $j$ considered in Figure 2, it can clearly be observed that not a single extreme outlier (one very small p-value corresponding to one individual subject) is responsible for the statistical significance with respect to the combination test statistic $T_j$ but that the combined information contributed by all $n$ subjects supports our statistical conclusions. Regional p-values $p_{\mathit{ij}}^{\mathit{\text{REGION}}}\equiv 1$ can occur, if none of the voxels belonging to region $j$ has been selected in Steps 1 and 2 described before for a certain subject $i$ . By construction of $T_j$ this essentially means that the “effective sample size” for such a region is reduced, while the number of degrees of freedom for the null distribution of $T_j$ remains unchanged.

5.1 Statistical sensitivity

This paper introduces a new strategy to analyze fMRI data and demonstrates its performance by drawing a comparison to two studies by Siegmund et al. In their second study, they investigated a new research question which was based on a design with reduced statistical power regarding the network of brain areas activated during bottom-up program comprehension. Only by using prior knowledge of the location of identified clusters from the first study and conducting a regions-of-interest analysis with increased statistical power, they were able to exceed standard thresholds of statistical significance but restricted to the subset of areas identified in the first study.

Research that similarly aims to detect small differences between cognitive processes with state-of-the-art methods would need to increase statistical power, e.g. by conducting two studies: first, to identify the network of brain areas involved in the overarching cognitive process, and second, to identify differential effects within the activated brain areas. The method of utilizing anatomically predefined regions of interest as determined in each individual brain is a candidate for a more sensitive analysis as compared to statistical testing of single voxels with identical coordinates of a template brain. This is possibly due to allowing more flexibility with respect to the exact anatomical location of functional units across different brains. Our example of using ROIs from the Harvard-Oxford atlas can easily be replaced in future applications by brain parcellation schemes that are based on more refined anatomical and/or functional databases. We demonstrated by simulations that our method can accurately detect signals based on region aggregation and outperforms standard mass-univariate approach while achieving a similar performance as cluster-based inferences without the need to define thresholds based on data smoothness via simulations. Further, we demonstrated that our method can find significantly activated brain areas without relying on prior knowledge on a real fMRI dataset. Moreover, additional brain areas were identified which have been described in two fMRI studies of programmers and interpreted as being involved in visuo-spatial processing. One study investigated manipulating data structures, which share similarities to spatial rotation (see Ref. 25) and one writing program code (see Ref. 27). Since the study by Siegmund et al. (2014) was the first study on program comprehension they used a rather small FDR corrected significance level ( $p<0.01$ as compared to the more common $p<0.05$ ). Possibly this is one reason why the two areas identified by the current approach were not identified there and as a consequence could not emerge in Siegmund et al. (2017). Thus, our current approach seems valuable in cases where new research questions are explored and little to no prior work is available and which would require careful statistical hypothesis testing to initially minimize false positives.

However, inferior frontal gyrus (BA 44) and inferior parietal lobule (with BA 40), shown to be significant in Siegmund et al., were not significant with our method. Therefore, we exemplarily investigated this difference with the activation cluster in BA 40. In Figure 5, we display a histogram of p-values for the (statistically significant) Aparc label in the middle temporal gyrus, in which a majority of the participants contribute with small p-values. In contrast, in Figure 6, we display a histogram of p-values for the inferior parietal lobule. Across the entire group, we also observe an accumulation of very small p-values, but only three of eleven participants contribute to this. Unlike the regions-of-interest analysis done by Siegmund et al., our method relying on combining Aparc p-values by Fisher’s method does not reject hypotheses with such a p-value distribution.

Figure 5. Histogram of p values of ctx_lh_G_temporal_middle label that is evaluated as significantly activated.

The plot without a heading is across all participants, while the eleven plots with headings show the distribution across individual participants. The majority of participants is contributing to the significance.

Figure 6. Histogram of p values of ctx_lh_G_pariet_inf-Angular that is not evaluated as significantly activated.

The plot without the heading is across all participants, while the eleven plots with headings show the distribution across individual participants. While there are many small p- values, only a few participants contribute these small values.

By transforming the Aparc labels into the Talairach space, we observed that the voxels included in the activation cluster of BA 40 do not perfectly align with the Aparc label that should cover this brain area, i.e. the inferior angular gyrus of the parietal lobe. Table 4 shows that the overlap to the activation cluster in BA40 is less than 50% and that only around 5% of its voxels are assigned to an Aparc label at all. In contrast, Table 3 shows that at least 75% of the voxels of the BA 21 cluster are aligned to Aparc labels of which two are significant with our method. Thus, a small activation cluster within a large anatomical region may get lost with our approach. We further discuss this potential drawback of the used anatomical segmentation and possible remedies in Section 5.3.

From the methodological point of view, our main contribution consists in a novel way how to combine evidence: Instead of aggregating single voxel data over all participants by mapping them to a standard brain template, we define subject-specific regions and combine the evidence on the level of these regions by means of a combination function. This is a generic methodological approach which is not restricted to the specific study setup of Siegmund et al., but can be applied to essentially any fMRI study design, involving an arbitrary number of contrasts.

Furthermore, also the (final) combination step of our proposed approach is generic in the sense that instead of the Fisher combination function any other (appropriate) combination function for p-values may be used. Recently, there has been a renewed interest in p-value combination methods; see, e.g., Ref. 52 (with discussion), Ref. 48, and Ref. 49.

5.2 Comparison to related fMRI analysis methods

Under the multiple testing framework, testing of grouped null hypotheses with (potential) application in fMRI is an active research topic. In Ref. 23 as well as in Ref. 4, clustering techniques were employed to define regions of interest, and the authors incorporated the heterogeneous cluster sizes in a weighting scheme for the linear step-up test from Ref. 6. In the same vein, the authors of Ref. 24 as well as Ref. 56 made use of the different proportions of true null hypotheses in each of the groups in their proposed weighting. A Bayesian variant of this idea has been derived in Ref. 30. Hierarchical methods, which exclude groups without strong evidence for the presence of signals in several stages of data analysis, have been worked out by several researchers; see, e. g., Refs. 3, 40 and 55. However, these methods rely on combining the subject-specific data on the voxel level, which is a standard technique as mentioned, for instance, in Section 5 of Ref. 29 Also on the basis of (combined) voxel data, a hierarchical independent component analysis for the comparison of brain functional networks has been proposed in Ref. 41. To the best of our knowledge, our proposed method to define subject-specific regions by means of regional labels and to combine the resulting subject-specific regional p-values by means of a combination function is a novel idea.

It is worth mentioning that both our proposed methodology and the methods from previous literature that we are considering in our comparisons make certain assumptions about the data: First, a linear relationship between the features encoded in the design matrix and the response is assumed. Thus, non-linear effects cannot (fully) be captured by models of the type (1). Second, Gaussianity of error terms is assumed, implying symmetry and light tails, among other things. It can be of interest to validate or to test these assumptions for a concrete dataset at hand. However, since this aspect is not the main focus of our present work, we defer such investigations to future research. In particular, modeling the data that we have re-analyzed in this work with different model assumptions would diminish the comparability of our data analysis results with the results from previous analyses. One simple (numerical) robustness analysis regarding the Gaussianity assumption can be performed by simulating data with a different error distribution and comparing the results with the results presented here. To this end, the source code which is available as supplementary material for this article may be helpful for practical implementation.

5.3 Outlook: from anatomical to functional aggregation

Instead of utilizing a single voxel GLM group analysis in common brain templates, we used a parcellation of the brain for each individual participant into regional labels, here Aparc labels, before aggregation into the group analysis. This procedure provides more labels than a traditional Brodmann atlas. Still, some of the regions are very large and thus presumably contain several functional areas. There is ongoing research to subdivide the brain based on cytoarchitectonic details, e.g. the Jülich-Brain (see Ref. 1). This will provide more detailed parcellation schemes in the future and that could easily be integrated into our proposed method. Another refinement could be to use functionally defined brain regions for our presented methodology. A study that implements functional localizers could identify participant-specific functional maps of the brain for well-defined standardized tasks (e.g., see Ref. 33). Then, in the analysis, our presented methodology can aggregate across all participant-specific brains with less imprecision than traditional methods. We would like to note that such functional localizers are currently restricted to research areas that include brain regions with specific functional specialization; cf. Refs. 20 and 39. The studies we presented in this paper are concerned with a rather complex cognitive task. Moreover, understanding program code is highly individual since different programmers rely on different comprehension strategies based on their preferences, domain knowledge, and experience (e.g., Refs. 8, 9 and 44).

Ethical considerations

Our paper did not acquire any new data with human participants. The original studies of Siegmund et al. were approved by the responsible ethics board of the University of Magdeburg (Application: 87/14).