|Home | About | Journals | Submit | Contact Us | Français|
When both structural magnetic resonance imaging (sMRI) and functional MRI (fMRI) data are collected they are typically analyzed separately and the joint information is not examined. Techniques that examine joint information can help to find hidden traits in complex disorders such as schizophrenia. The brain is vastly interconnected, and local brain morphology may influence functional activity at distant regions. In this paper we introduce three methods to identify inter-correlations among sMRI and fMRI voxels within the whole brain. We apply these methods to examine sMRI gray matter data and fMRI data derived from an auditory sensorimotor task from a large study of schizophrenia. In Method 1 the sMRI–fMRI cross-correlation matrix is reduced to a histogram and results show that healthy controls (HC) have stronger correlations than do patients with schizophrenia (SZ). In Method 2 the spatial information of sMRI–fMRI correlations is retained. Structural regions in the cerebellum and frontal regions show more positive and more negative correlations, respectively, with functional regions in HC than in SZ. In Method 3 significant sMRI–fMRI inter-regional links are detected, with regions in the cerebellum showing more significant positive correlations with functional regions in HC relative to SZ. Results from all three methods indicate that the linkage between gray matter and functional activation is stronger in HC than SZ. The methods introduced can be easily extended to comprehensively correlate large data sets.
Developments in neuroimaging techniques, especially in magnetic resonance imaging (MRI), have enabled researchers in the past decade to perform in vivo studies of the structural and functional brain. Structural MRI (sMRI) is widely employed to image different tissue distributions and functional MRI (fMRI) is used to identify brain activation regions for a certain task. During an fMRI experiment a subject is asked to perform a task while the scanner records the blood oxygen level dependent (BOLD) changes of different brain regions. fMRI has become a popular functional neuroimaging technique since it is less invasive, less expensive and has better spatial and temporal resolutions than previous techniques such as positron emission tomography (PET). sMRI and fMRI allow investigators to study groups of patients and healthy controls and identify differential brain characteristics between groups. The identified anatomical and functional measures or a connection between them may be subsequently used as indicators of susceptibility to the disorder.
In most fMRI studies sMRI is also acquired. These two modalities provide unique information about the brain; however, their analysis is typically performed separately. The sMRI images are often used to co-register different subjects’ brains onto a common template (Toga and Thompson, 2001) or used as a rendering surface to visualize overlaid functional activation (Corbetta et al., 1998). The brain is a vastly interconnected organ and it is reasonable to expect that changes in local morphological structures may result in modulations of brain activity in distant regions (Mesulam, 1998). The correlation or the connection between structural regions and functional activation of the whole brain has not been well established. Studies have been designed to examine localized correlations between gray matter volumes with fMRI or to correlate structural and functional data from certain regions in the brain (Hasnain et al., 1998, 2001; Siegle et al., 2003). Existing tools for examining joint information include region-based approaches such as structural equation modeling (McIntosh and Gonzalez-Lima, 1994) and dynamic causal modeling (Friston et al., 2003). The above methods were primarily designed to investigate functional connectivity between brain regions. These studies and techniques are important and needed; however, they are limited by a priori knowledge of the implicated brain regions and lack detailed examination of the relationship between distributed brain voxels. The confound of anatomical differences or registration error contributing to differential functional activation was addressed by integrating differences in underlying gray matter as a nuisance variable into the general linear model (Oakes et al., 2007). This method addresses an important issue but investigates functional activation differences due to underlying structure only within the same voxel and hence does not examine relationships between structure and function at distant brain regions. There are other methods (Fan et al., 2007; Ford et al., 2002) that use structural and functional data for classification purposes; however, they do not examine the inter-relationship between structurally and functionally derived data.
In our approach we assess cross-correlations among all sMRI and fMRI voxels. This task is straightforward for a small numbers of voxels, but becomes challenging for whole brain studies due to the magnitude of necessary computations and temporary storage needed. With tens of thousands of voxels from two modalities to be comprehensively correlated a cross-correlation matrix with billions of values needs to be constructed. In this paper we introduce techniques to iteratively compute the large cross-correlation matrix and reduce it to statistical measures that can be used to identify discriminating features between healthy and patients groups. We have used a related approach to evaluate function–function linkages (Michael et al., in press). The correlation for each pair of sMRI–fMRI voxels is computed across subjects belonging to a certain group. Our method is a data driven approach with no prior assumptions about the brain regions to be investigated.
We apply our methods to a group of patients with schizophrenia (SZ) and healthy controls (HC). Schizophrenia is a complex and chronic mental disorder that affects about 1% of the population across all cultures, sexes and socioeconomic groups (Bhurga, 2005). The clinical heterogeneity of schizophrenia, characterized by a myriad of symptoms and signs, presents a difficult target for research into its neurobiological mechanisms (Andreasen, 2000). No neurobiological marker of the disorder is present and no specific brain region has been found to be solely responsible for the disorder. Schizophrenia is currently diagnosed on the basis of patient’s self reported experiences and observed behavior.
Studies separately analyzing structural and functional images have found that multiple brain regions appear to be affected in schizophrenia (Goldstein et al., 1999; Honea et al., 2005; Niznikiewicz et al., 2003). Structural findings suggest enlarged ventricles, reduced volume of temporal lobe and superior temporal gyrus, with moderate evidence for frontal lobe volume reduction, and some evidence for persistent cavum septi pellucidi and abnormalities in basal ganglia, corpus callosum, thalamus and cerebellum. Functional results suggest abnormal connectivity between temporal and frontal brain regions and aberrant activation in the dorsolateral prefrontal cortex (DLPFC) (Manoach et al., 2000; Weinberger et al., 1986). An overwhelming number of studies report that for a given task, abnormal function is not localized to a specific brain region, but instead the abnormalities involve a network of brain regions.
Schizophrenia is likely associated with disruptions in the connectivity between different cortical regions (Andreasen et al., 1999;Andreasen et al., 1998; Breakspear et al., 2003; Friston, 1998; Kubicki et al., 2007; Lim et al., 1999; Ross and Pearlson, 1996). Studies that investigate the linkage between structure and function while taking into account all brain voxels are thus needed to investigate the differences between HC and SZ. The term ‘linkage’ in this paper is used to denote the relationship, association or correlation between two variables. A joint analysis of this nature enables one to do a comprehensive study of the interactions between sMRI and fMRI. It is reasonable to expect that large functional networks are associated with features of gray matter but to different degrees in HC and SZ. The disconnection hypothesis of schizophrenia (Friston, 1998) states that neural mechanisms in schizophrenia are not localized to any one area alone, but it is the integrity of the interconnections between brain regions that is compromised. Thus, we predicted that the relationship between structure and function will be weaker in SZ than in HC.
sMRI data were derived from T1-weighted scans and with voxel-based morphometry (Ashburner and Friston, 2000) the images were segmented to obtain gray matter (GM) concentration. GM concentration is the percentage of GM content at a voxel location. By form, given in the title, we mean the tissue composition of a voxel. We selected GM, since GM brain regions are the information processing centers of the brain. Many previous studies have found gray matter volume and density differences in SZ and motivated us to investigate its correlation with function (Ananth et al., 2002; Gur et al., 1999;Hulshoff Pol et al., 2001). fMRI data were acquired while subjects performed an auditory sensorimotor (SM) task and were reduced to activation maps by regressing the task time course and other models against the data using the GLM. The SM task was designed to activate the auditory cortex robustly since previous studies (Calhoun et al., 2007; Kiehl and Liddle, 2001; Machado et al., 2007; Schroder et al., 1999) have shown abnormal activation in auditory temporal lobe regions in SZ.
Subjects recruited for this study were scanned at four different sites: the University of Iowa Hospital (IA), Harvard’s Massachusetts General Hospital (MA), the University of Minnesota (MN) and the Mind Research Network (NM). Patients were recruited from inpatient and outpatient psychiatric clinics, group homes, referrals from physicians, advertisements and by word-of-mouth. Controls were recruited from advertisements, fliers and word-of-mouth. All participants provided written, informed, IRB approved consent at their respective sites and were compensated for their participation. Subjects analyzed in this study had normal hearing (assessed by self report) and were able to carry out the fMRI task. A total of 70 healthy controls (HC) and 70 patients with chronic schizophrenia (SZ) are used in this analysis. We checked that the structural brains and functional activation regions were consistent across subjects through a correlation analysis. The two subject groups are matched for age, sex and parental socioeconomic status and their demographics are presented in Table 1. In Table 1 all demographics except years of education have insignificant group differences. A large number of subjects have to be removed to match all between group subject demographics. To avoid this we select a subset of subjects to match education between the two groups and repeat the tests to show the robustness of our result. The significance of structural–functional correlation between the two groups is comparable since the number of subjects in each group is equal.
HC in the study were screened to ensure they were free from diagnostic and statistical manual of mental disorders (DSM-IV) axis I or axis II psychopathology, assessed using a modified version of the comprehension assessment of symptoms and history (CASH) (Andreasen et al., 1992). They were interviewed to determine that there was no history of psychosis in any of their first-degree relatives. Patients met criteria for schizophrenia based on structured clinical interview for DSM (SCID) or CASH and were confirmed by review of their case file. Average patient symptom measures, positive (Andreasen, 1984), negative (Andreasen, 1981) and disorganization, are reported in Table 1. Out of the 70 SZ, 58 were receiving the following antipsychotic medication (number of subjects in parenthesis): clozapine (12), quetiapine (8), olanzapine (9), ziprasidone (4), risperidone (9), aripiprazole (9), quetiapine and aripiprazole (1), perphenazine (1), haloperidol (3), thiothixene (1), fluphenazine and thiothixene (1). Information was not available for 12 SZ.
fMRI data were collected while subjects performed an auditory sensorimotor (SM) task. Prior to the scan all subjects were instructed until they were able to perform the task adequately. The participants were requested to actively participate and respond as quickly and accurately as possible. The SM task was designed to robustly stimulate the auditory cortex. Sound stimuli were incorporated into E-prime scripts (http://www.pstnet.com) run on a Windows machine and were presented via sound insulated, MR-compatible earphones (Avotec, Stuart, FL). The task consisted of an on/off block design (Fig. 1). During the on-block, 200 ms tones were presented with a 500 ms stimulus onset asynchrony (SOA). There were eight tones at different pitches along a scale and they were presented in ascending and descending cycles for the duration of the ‘on’ block. No tones were presented during the ‘off’ block. Upon each presentation of the tone a subject was required to make a speeded right thumb button-press response through an MR-compatible input device (https://www.mrn.org/facilities/mind-input-device). There were two runs, each lasting about 4 min.
A T1-weighted sMRI was acquired at each site. For MA and NM the scans were acquired on a Siemens scanner at 1.5 Tesla (T) with the following parameters: repetition time (TR)/echo time (TE) = 12 ms/ 4.76 ms, bandwidth = ±150 kHz, flip angle (FA) = 20°, slice thickness= 1.5 mm, voxel size = 0.63 × 0.63 × 1.5 mm3, field of view (FOV)= 161 mm, and pulse sequence = gradient echo. For IA the scans were acquired on a GE Signa scanner at 1.5 T with the following parameters: TR/TE = 20 ms/6 ms, bandwidth = ±180 kHz, FA = 30°, slice thickness= 1.6 mm, voxel size = 0.66 × 0.66 × 1.6 mm3, FOV = 168 mm, and pulse sequence = gradient echo. The MN scans were acquired on a Siemens Trio at 3 T with the following parameters: TR/TE = 2530 ms/3.8 ms, bandwidth = ±180 kHz, FA = 7°, slice thickness = 1.5 mm, voxel size = 0.625 × 0.625 × 1.5 mm3, FOV = 160 mm, and the pulse sequence = MP-RAGE. At all four sites brain images were acquired in a coronal orientation.
Functional data were acquired at all four sites with EPI sequences on Siemens 3.0 T scanners, except at the NM site where a 1.5 T Siemens scanner was used. Data were collected from each participant while performing the SM task. The parameters for the functional scan are as follows: TR/TE = 2 s/30 ms (TE = 39 ms for NM), bandwidth (BW) = ±100 kHz= 3126 Hz/pixel, FA = 90°, slice thickness= 4 mm, gap between slices=1 mm, voxel size=3.4×3.4×4 mm3, FOV= 22 cm, pulse sequence = PACE-enabled, single shot, single-echo echo planar imaging (EPI), scan plane = oblique axial, AC-PC; acquisition matrix=64×64, number of slices=27, ascending sequential acquisition.
For MA and NM, three T1’s were coregistered to each other and an average T1 was computed for segmentation and smoothing. IA and MN acquired only a single T1 image. Tissue classification, bias correction, image registration and spatial normalization were performed using voxel-based morphometry (VBM) (Ashburner and Friston, 2005) in SPM5 (http://www.fil.ion.ucl.ac.uk/spm/software/spm5) where the above steps are integrated into one model. Unmodulated normalized parameters were used for segmentation as previously applied in two large VBM studies (Meda et al., 2008; Segall et al., 2008) to segment the brain into white matter (WM), gray matter (GM) and cerebral spinal fluid (CSF) probabilistic maps. For each brain voxel the three concentrations from these maps add up to one. After segmentation, the GM images were smoothed with a Gaussian kernel of full-width half maximum (FWHM) of 10 × 10 × 10 mm3. The voxel size of all images was resliced to 3×3×3 mm3.
The SPM5 software package was employed to perform fMRI preprocessing. Slice timing was performed with the middle slice as the reference frame. Images were realigned using INRIalign, a motion correction algorithm that is unbiased by local signal changes (Freire and Mangin, 2001; Freire et al., 2002). Data were then spatially normalized into the standard Montreal Neurological Institute (MNI) space (Friston et al., 1995) with affine transformation followed by a non-linear approach with 4 × 5 × 4 basis functions. Images (originally collected at 3.4 × 3.4 × 4 mm3) were then slightly downsampled to 3×3×3 mm3, resulting in a data cube of 53 × 63 × 46 voxels. Finally data were spatially smoothed with a Gaussian kernel of FWHM of 9×9×9 mm3 (White et al., 2001).
The SM task was modeled as a block design as shown in Fig. 1. The task regressor was created by modeling the on–off blocks as rectangle functions convolved with the default SPM5 canonical hemodynamic response function (HRF). Drift was modeled by a high pass filter with cutoff at 128 s. The contrast map was computed with these regressors using GLM in SPM5. Contrast maps from the two runs of the task were averaged for each subject to create the final map input into our approach.
The brain maps contained about 153 k (53 × 63 × 46) voxels and after the exclusion of non-brain voxels the final number of voxels was about 60 k (N). The reduced three dimensional voxels were vectorized along different columns of a matrix where subjects belonging to a certain group (SZ or HC) were stacked along the rows (Fig. 2). Let Xi and Yj be the column vectors across all subjects for the ith voxel from sMRI (gray matter concentration map from VBM) and the jth voxel from fMRI (activation map from GLM), respectively. Our interest was in finding the correlation (using Eq. 1) between Xi and Yj where i and j varied independently from 1 to N.
A cross-correlation matrix (RSF) of size N × N is required to store all structural–functional correlations and making such a matrix is not easy due to limitations in computer memory. Even if RSF is computed its interpretation will be based on some reduced statistics. We introduce three methods to find statistics by iteratively computing RSF.
The histogram of all the elements of RSF gives a general idea of how the structural–functional correlations are distributed. Steps to find the histogram of RSF are listed below (see Fig. 3).
This method essentially finds the histogram of RSF (N × N) by iteratively finding the histograms of its rows and finally adding the histograms. In the next paragraph we show that the histogram computed through this method and the histogram of all the elements of RSF is identical given that the binning intervals are identical. In that sense we do not lose information but we are reducing N × N correlations of the RSF matrix to 100 numbers of occurrences, a form that can be easily represented and analyzed.
In Fig. 4 the histograms of each row of RSF are shown as number of occurrences within different binning intervals. The lower (L1), upper (L2) bounds and the bin interval (Δx) are kept constant across the rows of RSF. In our application L1 and L2 are equal to −1 and +1, respectively, since they are the lower and upper bounds of the correlation coefficient and Δx equals (L2–L1)/100. The number 100 was selected arbitrarily. Hi is the histogram of the elements of the ith row of RSF and fij is the number of occurrences in the jth bin of Hi. Histogram Hi is essentially a collection of occurrences *(fi1, fi2, …, fij, …, fi100). If L1, L2 and Δx are kept constant and if RSF is fully computed the number of occurrences within a certain bin will be the summation of occurrences of the rows of RSF within the same bin. Hence the summation of Hi is identical to the histogram of RSF if it had been fully computed.
In Method 1 we find the distribution of structural–functional correlations but the spatial information of where these correlations are located is lost. In Method 2 we try to retain average spatial distribution of the correlations. This is achieved by collapsing either the rows or the columns of RSF (Fig. 3b) and then reconstructing the p-values on to a brain map.
If RSF is collapsed along the columns, the p-value at a certain voxel in the reconstructed brain map corresponds to the significance of the mean correlation of the structural voxel at that location with all functional voxels being non-zero. If collapsed along the rows the p-value at a certain voxel corresponds to the significance of the mean correlation of the functional voxel with all functional voxels being non-zero. To compute the p-value we performed a one sample t-test and computed the t-value using the mean and standard deviation of the correlations at that voxel. The number of subjects (=70) was used as the original degrees of freedom (dof) since the correlations were computed from 70 independent subjects. To examine how each structural voxel is correlated to all functional voxels a total of N maps are needed and inspecting such a large number of brain maps is not an easy task. It should be noted that this method gives an average sense of the correlations. In this method we are reducing the N × N correlations of RSF to a map of N p-values of correlations.
In method 2 we retained average spatial correlations but did not retain information on structural–functional inter-regional correlations. With Method 3 it is possible to find structural–functional regions that are significantly correlated.
In Fig. 3c, various segments are represented by different colors. The orange segment, for example, corresponds to how structural voxels from the 1st atlas region correlate to functional voxels from the 2nd atlas region. In this method we are reducing the N × N correlations of RSF to 116× 116 inter-regional p-values of correlations. Here the p-value indicates the significance of the mean correlation of a segment being non-zero.
Method 1 was applied to the group of HC and SZ separately and the structural–functional correlation histogram was computed with all brain voxels. In Fig. 5a the y-axis represents the number of occurrence of different correlation values given along the x-axis. The summation of the number of correlation occurrences in the histogram equals the number of elements of RSF (N × N). A mean correlation value of zero in the histogram indicates that on average the structural gray matter concentration (across subjects) and functional activation values (across subjects) are uncorrelated. A positive correlation means that an increase or decrease in gray matter concentration corresponds to an increase or decrease in functional activation respectively. A negative correlation value conveys that a change in gray matter concentration has a change in functional activation in the opposite direction.
In Fig. 5a the SZ histogram is shown in red and the HC in blue. The difference between SZ and HC histograms is not very clear when they are drawn on the same plot. To enhance the difference in Fig. 5b the difference histogram (HC–SZ) is illustrated. From this plot it is evident that the number of zero correlations is higher in SZ than in HC and the number of +0.2 and −0.2 correlations is higher in HC than in SZ. To find the significance of the correlation between functional activation and gray matter concentration we converted the correlation value (±0.2) that showed group difference into a t value using Eq. 2, where r = ±0.2 and ndof = 70 (number of subjects). We then converted the t value to its corresponding p-value. After these steps for an r = ±0.2 we obtained a p-value of 0.048. This value does not indicate the significance of the difference between the two groups but the significance of correlation (at r = ±0.2) between structural and 436 functional data within a group is p<0.05.
The result presented above was based on a cohort of subjects with insignificant group differences in age, sex and parental socioeconomic status but a significant group difference in years of education. To match education we selected a subset of 65 subjects from each group. The mean, standard deviation and the significance of difference between the two groups are as follows: HC = 14.6 ± 1.5, SZ = 14.1 ± 2.4, t-value = 1.4 and p-value = 0.15. The structural–functional correlation histograms (as in Fig. 5) were computed and the result that HC having higher correlations than SZ was consistent with this subject group as well.
In the previous analysis structural–functional correlations were found using all brain voxels. During the preprocessing steps we had removed voxels that were ‘noisy’ as a result of scanner issues. Another set of ’noisy’ voxels appear when voxels do not activate consistently across all subjects. We were interested in checking the above result for functional voxels with significant activation and structural voxels with high gray matter concentration. Our test was to find if the pattern of HC having higher correlations than SZ changed across different levels of significant voxels. Functional voxels that had the highest group (HC + SZ) mean t-values for the SM task and structural voxels with the highest group mean gray matter concentration were used for this analysis. The voxels that pass a certain threshold in the group mean map were chosen from each subject to compute the histogram. We compute the histograms for different combinations of structural and functional thresholds and the difference histograms (HC–SZ) are presented in the rows of Fig. 6a. In Fig. 6a the difference histograms (HC–SZ) are shown as a colormap. Shades of red correspond to HC having higher number of correlations than SZ and shades of blue correspond to SZ having higher number of correlations than HC at correlation values given by the x-axis. In Fig. 6a SM activation map t-value is changed (between 0 and 4.0 in steps of 0.25, separated by dashed lines) in the outer loop and GM concentration is changed (between 0 and 0.8 at steps of 0.1) in the inner loop. The first row in Fig. 6a corresponds to all structural and functional voxels and the last row to structural voxels with GM concentration above 0.8 and functional voxels with t-values above 4.0
At lower SM t-value thresholds (0.0–1.5) and lower GM thresholds (0.0–0.5) HC show higher number of correlations around +0.2 and –0.2 and SZ show higher number of correlations around zero. This pattern is consistent with the histogram found in the previous section using all brain voxels. At higher GM thresholds (0.5–0.8), independent of SM t-values, HC show higher number of positive correlations and SZ a higher number of negative correlations. This is seen in Fig. 6a, as GM thresholds are increased from 0.5 to 0.8, blue and red shades move to the left and right of correlation =0 axis, respectively. In Fig. 6a color intensity reduces at higher GM and SM thresholds. This is due to the fact that at higher thresholds a lower number of voxels are selected to compute the histogram and hence the number of occurrences is also lower. To find the significance of the histograms at different thresholds we normalize the histogram by dividing it by the total number of correlation occurrences. A normalized histogram is essentially a probability density function (pdf). In Fig. 6b the difference between HC and SZ normalized histograms are illustrated in a manner similar to Fig. 6a. From Fig. 6b we see that a more significant difference between the two groups occurs at higher SM and GM thresholds. At these thresholds controls show more positive and patients more negative structural–functional correlations.
In this section we present results, using Method 2, where the spatial content of the correlations is retained. In Fig. 7 we present the p-values of how a functional voxel, on average, is correlated to all structural voxels. This p-value is presented at the brain location of the functional voxel. In Fig. 7a and b the p-value indicates the significance of the mean correlation of a functional voxel with all structural voxels being non-zero. The p-values in Fig. 7c correspond to the significance of the mean correlation of one group being higher than the other. Voxels with p-values less than 0.01 are mapped in Figs. 7a–c for HC, SZ and HC–SZ, respectively. In Figs. 7a and b if the mean correlation is positive the p-value is denoted with a shade of red or yellow and if the mean value is negative it is denoted by a shade of blue. All reported p-value regions are corrected for multiple comparisons using Bonferroni correction.
In HC there were large regions of positive correlations in the following anatomical regions: superior temporal gyrus, pre- and postcentral gyri, superior frontal gyrus and insula. Negative correlations were seen in the following regions in HC: middle temporal gyrus, inferior temporal gyrus and middle occipital gyrus and cuneus. In SZ the following regions showed positive correlations: inferior parietal lobule, lentiform nucleus, middle frontal gyrus, superior temporal gyrus and insula. Negative correlations in SZ were seen in precuneus, middle occipital gyrus, cuneus, fusiform gyrus, cingulate gyrus, middle temporal gyrus, culmen and anterior cingulate. In Table 2 we present regions that show significant difference (p<0.01) between HC and SZ of how a certain functional voxel, on average, is correlated to all structural voxels. These regions were converted from MNI space to Talairach coordinates and entered into a database (http://ric.uthscsa.edu/projects/tdc/) to provide the labels and volumes of contiguous regions.
InFig. 8 we present the p-values of how a structural voxel, on average, is correlated to all functional voxels. These p-values were found in the same manner as explained for the previous result. In Figs. 8a–c regions with p<0.01 (after Bonferroni correction) are presented for HC, SZ and HC–SZ, respectively.
An interesting result was that in HC, cerebellar regions showed more positive correlations and prefrontal regions showed more negative correlations than did SZ. From Figs. 8a and b it is seen that HC show larger and more significantly correlated regions than SZ. The following brain regions showed positive correlations in HC: cerebellum, middle occipital gyrus, pyramis, superior temporal gyrus, inferior semi-lunar lobule, precuneus, thalamus and cuneus. In the following brain regions, negative correlations were seen in HC: medial frontal gyrus, superior frontal gyrus, middle frontal gyrus, precentral gyrus, inferior frontal gyrus, cingulate gyrus, anterior cingulate and insula. SZ did not exhibit large contiguous regions of positive nor negative correlations. In Table 3 we present the regions and their volumes, that show significant (p<0.01) group differences of how a certain structural voxel is correlated to all functional voxels.
In this section using Method 3 we identify regions in the brain that have significant structural–functional linkages. With this method the structural–functional cross-correlation matrix (RSF) is reduced to a matrix of p-values of size 116 × 116 indicating the significance of nonzero mean inter-regional correlations. The cerebellar vermis from the structural data had significant (p<0.01, after Bonferroni correction) positive correlation in HC with functional data from the following regions: calcarine, cuneus, lingual gyrus, paracentral lobule and Heschl’s gyrus. In SZ there were no positive correlation below a p-value of 0.01 (after Bonferroni correction) but the basal ganglia from the structural data had significant negative correlation with the posterior cingulate from the functional data.
Our next interest was in determining inter-regional correlations that were significantly different between HC and SZ. We used two sample t-test to compute the t-value corresponding to a certain segment and then converted it to a p-value. Structural regions in the cerebellum had significant (p<0.01, after Bonferroni correction) positive correlation (HC>SZ) with functional regions in the calcarine, lingual gyri, occipital lobe and Heschl’s gyrus.
The main purpose of this paper is to present methods to fuse structural and functional data through spatial sMRI–fMRI cross-correlations across the whole brain. These methods enable one to compute statistics of the very large sMRI–fMRI cross-correlation matrix (RSF) through iterative computing. Our interest was to investigate how local structure (in our application gray matter concentration) may influence functional activity (probed through a sensory motor task) across the whole brain. This objective is motivated by the fact that the brain is a highly interconnected organ and also to check if the relationship between structure and function is compromised in healthy versus diseased brain (in our application schizophrenia).
The structural–functional connection is tested by finding the correlation between two vectors. The first vector is gray matter concentration of a brain voxel and the second functional activation from another brain voxel. These vectors are obtained across subjects belonging to a certain group (controls or patients). A strong correlation between these two vectors indicates that structural and functional variations across subjects are related. We repeat this correlation process for all structural and functional voxels to compute N × N structural–functional correlations, where N is the number of brain voxels. For easier analysis and better representation, the N × N correlations are reduced to metrics through three different methods we introduce with increasing levels of anatomical specificity. In Method 1 we obtain a histogram of the N × N correlations by reducing them to 100 numbers of occurrences at correlation values between +1 and −1. A higher number of zero correlations in one group indicate that the structural–functional correlation of the whole brain across subjects is weaker in that group. With this method we are ableto demonstrate the distribution of the correlations; however, the spatial location of the correlations is not retained. In Method 2 we preserve the average spatial maps of how voxels from Modality1 are correlated to all voxels from Modality2. By collapsing across the columns of the cross-correlation matrix we find how a certain structural voxel is correlated to all functional voxels and the corollary is obtained by collapsing across the rows. If a structural voxel shows high mean correlation with all functional voxels it indicates that the gray matter concentration variation of that voxel across subjects and functional activation of all brain voxels across subjects on average are highly related. With Method 2 the average correlation is found across the whole brain and useful information can be lost. For example some regions in the brain can have negative correlations and others positive correlations and an average can indicate a correlation value close to zero. This issue is partially addressed by reporting the significance of non-zero correlation through p-values that accounts for this variation. We also report the significance of positive and negative correlations separately as different colors in Figs. Figs.77 and and8.8. In Method 3 we made an effort to find significantly correlated links between structural and functional brain regions as defined by the AAL atlas. The methods introduced are simple ways to reduce a large data set and answer three different questions about the large cross-correlation matrix: distribution of correlations, average spatial maps of correlations and average inter-regional correlations. The methods can be further improved to answer different questions. For example which pair of voxels has the most significant correlations, which segments of RSF have the highest positive or negative correlations or how two regions of interest are correlated. Depending on the nature of the question, reduction steps need to be devised and methods to represent the results also need to be developed. What we demonstrate in this work is that even using simple reduction techniques, it is possible to find features that significantly differ between groups of patients diagnosed with schizophrenia and healthy controls.
The approaches were applied to a group of patients with chronic schizophrenia and matched healthy controls to investigate if the degree of association between the gray matter concentration and functional activation probed through an auditory sensorimotor task is different in these two groups. Our results reveal several interesting findings. The main finding was that patients with schizophrenia show less correlation across subjects between structural data (gray matter concentration) and functional data (obtained from the activation map of a sensorimotor task) than healthy controls. This result was confirmed with all three different methods and was found to be consistent. The results reported in this study were derived from a total of 140 subjects (70 in each group). The large number of subjects is an advantage as it increases the statistical power of the analysis.
In Method 1 we found that at correlation values centered around +0.2 and −0.2 the controls have a higher number of correlations than patients and at zero correlation patients had a higher number of correlations than controls. This finding suggests that the connection between gray matter concentration and functional activity as probed by the sensorimotor task is weaker in patients than in controls.
Method 2 results showed contiguous regions of high correlations in both patients and controls. This indicates that highly correlated regions were not located in randomly distributed brain regions, but were spatially clustered, and increases the validity of our results. In the map of how a functional voxel, on average, is correlated to all structural voxels (Fig. 7) both patients and controls show positive correlations in the superior temporal gyrus but healthy controls show larger and more significant correlation regions than in patients. This result indicates that the variation of functional activity in the superior temporal gyrus across subjects and variation of gray matter concentration of the whole brain across subjects was more correlated in controls than in patients. The question of whether this difference was caused by inaccurate registration of subject data is reasonable. The spatial normalization step during the preprocessing step used the same template to normalize patient and control subject data. Minimal topological differences can be attributed to spatial normalization but volume differences as large as listed in Table 2 are unlikely to be due to registration issues. Our results indicate functional regions that were not directly involved in the functional task to be having high correlations with structural regions. This result was made possible since our methods did not make prior assumptions on regions of interest (ROI). In Fig. 8 and Table 3 controls demonstrate significantly more correlated regions than patients. The main findings included cerebellar regions in controls had positive correlations and prefrontal regions had negative correlations. In controls the correlations in cerebellar regions in the left hemisphere were larger and stronger than the corresponding right region (see Fig. 8a) even though higher activations were observed in the right cerebellar regions (Fig. 1). We performed a two sample t-test on gray matter concentration but did not find cerebellar regions to be significantly different between the two groups. This indicates that between group differences in gray matter concentration may not have contributed to the differences in correlations.
In Method 3 we computed structural–functional inter-regional correlations and did not see significant correlations between the same structural and functional regions. This should encourage researchers to investigate correlations between distant brain regions. This point supports the hypotheses that local morphological structures can have correlations with distant functional activation.
A possible reason for differences in correlation can be attributed to differences in activation regions between the two groups. The activation maps for the functional task for the patient and control groups are presented in Fig. 1. The maps indicate that both patients and controls activated similar functional regions and hence did not possibly contribute to differences in structural–functional correlations. We also checked how our correlation result changed over different threshold for gray matter and functional activation. In Fig. 6 it is seen that the result of patients having higher number of zero or negative correlations is consistent across all combinations of GM and functional activation thresholds.
It is possible that higher correlations can be introduced as a result of head motion while correlating voxels within the same brain. Let us assume that during data acquisition group 1 has still heads inside the scanner and group 2 has head motion. If in group 2 a certain voxel is active, due to head motion, its activation will spread to neighboring voxels in the activation map. This can make correlations higher in group 2 than in group 1. In our results we see less correlation in patients than in controls. However, in reality it is the patients who have more head motion or tremors than the controls. If head motion had contributed to the results then our results should have been flipped with patients showing more correlation than controls.
In the theory of cognitive dysmetria (Andreasen et al., 1998) it is hypothesized that the system that is disturbed in schizophrenia is more distributed and it is modeled as a dysfunction in cortical–subcortical–cerebellar circuitry. This theory states that deficits in the complex distributed neural circuitry within the brain can express itself with a broad range of symptoms. The cerebellum has substantial connections with prefrontal cortex and can perform parallel processing supported by its array structure and large number of condensed cells (Andreasen et al., 1998). In our results we report regions that are widely distributed across the whole brain, especially in the cerebellar, subcortical and prefrontal regions, which interestingly are the nodes that fall under the rubric of the cognitive dysmetria model.
The results reported in this study also support the disconnection hypothesis (Friston, 1998) of schizophrenia. The disconnection hypothesis states that the disconnection is explicitly functional, not anatomical. This effective connectivity, as opposed to anatomical connectivity, is hypothesized as a result of differences in synaptic efficacy (Friston, 1998). In this work we did not directly investigate the anatomical connections between different regions of the brain but the correlation between gray matter concentration and functional activation. This investigation does not report on synaptic efficacy but provides evidence that further sources can exist for sMRI/fMRI linkage differences in schizophrenia.
There are several advantages of the methods we have introduced. The approaches are simple, easy to implement and efficient. They are simple, since they do not use complex mathematical equations, but rather the fundamental correlation equation to make all computations. They are easy to apply since the same simple computation is performed iteratively. In less than 40 min it was possible to compute the histogram for all possible combinations of all brain voxels (presented in Fig. 5) with MATLAB (R2006a) on a desktop machine (1.86 GHz CPU speed and 2 GB RAM). The methods are data driven and results obtained are after a comprehensive correlation analysis. We did not have prior hypotheses on which regions would show high correlations. Our results indicate non-task related regions (example left cerebellum) to be showing high correlation difference between the two groups. This result shows how data driven approaches can discover new features of the disorder. Another advantage is that the two modalities that are to be fused can have different resolutions or matrix sizes and thereby the loss of finer resolutions acquired can be avoided. In our application we matched the spatial resolution of sMRI and fMRI but this is not a necessary requirement. If the resolutions of the two modalities are different the cross-correlation matrix will be a non-square matrix. In data fusion the correct method to normalize the data sets to be fused is often unclear. In our methods we combine two data sets using correlation, a simple approach that also normalizes the two vectors. Correlation is also insensitive to the amplitude of the signal. For example, if there exists a relationship between non-task related functional activation and gray matter concentration, correlation will be able to pick it up. Another advantage of the method introduced is that it reduces issues pertaining to brain registration (Thompson et al., 2000). Instead of performing ROI analyses our method investigates voxels from the whole brain. In ROI methods perfect registration is crucial and even the slightest misalignments can cause errors in results. Further improvements can be made to the methods introduced to fuse data from other modalities, for example, to EEG or genetic data with sMRI or fMRI.
Further studies are needed to confirm our findings that the linkage between structure and function is weaker in patients with schizophrenia than in healthy controls. Our analysis was based on just one functional task and it will be important to confirm the results with other tasks that activate similar and different cognitive processes. It may be useful to apply this method on related and unrelated fMRI tasks to replicate and to assess the generalization of the findings. We measure brain function through fMRI, an imaging modality that is being widely used to measure brain function. fMRI captures the blood oxygen level dependent (BOLD) signal and is an indirect measurement of neural activity but has inherent spatial and temporal limitations (Menon and Kim, 1999). EEG and MEG are other modalities that can capture brain function but is not appropriate for this study due to their low spatial resolution. Our functional maps were obtained through GLM, a model based approach, where the hemodynamic response function is assumed to be identical for all subjects at all time points for all brain regions. There are reports that in schizophrenia there can be reduction in fMRI amplitude response (Carter et al., 2001), delay in HRF (Ford et al., 2005) and that HRF can vary from region to region (Miezin et al., 2000). Our method is insensitive to fMRI amplitude variations as explained above, but in this study we did not investigate the effects of HRF delay or altered HRF. However the between group differences we report are unlikely due to suboptimal HRF since the activation regions (Fig. 1) of the two groups are similar in overall activation. Independent component analysis (ICA) (Calhoun et al., 2001) is a model free method and can be used to find functional components that can vary as a result of changes in HRF. Subject demographics were not incorporated in our analyses; however, we show that our results were consistent when age, sex and education were matched between the two groups. Subject age had a range of about 40 years and it is possible that age-related differences exist. The groups included a mix of different handedness and medication, which were not included in the model. In order to assess specific demographics, a larger number of subjects would be required. In this work we also did not investigate how the structural–functional correlation would change with severity of symptom scores of schizophrenia. Our methods can be modified to check such a hypothesis. In our work we focused on the gray matter content in a brain voxel. Cortical thickness, gyrification or the physical connections between brain regions are other structural forms that may contribute to functional differences. In a previous study (Harris et al., 2007) it is found that prefrontal gyrification can predict those who develop schizophrenia. Prefrontal regions do show differences in our results as well, although, we did not investigate the relationship between gray matter concentration and gyrification or cortical thickness. In our method we attempt to identify linearly related linkages and in future work we intend to extend the methods to identify non-linear relationships, for example, using mutual information where higher order statistics can be investigated. A drawback of our method is that to find structural–functional relationships a group of subjects is needed. With our methods it is not possible to find the structural–functional relationships on an individual subject. This drawback makes it difficult to develop an algorithm that can classify a new subject as a patient or control. However, in a previous effort (Michael et al., 2008) we used the direction of shift of the functional– functional histogram (similar to Method 1) as a feature to classify when a new subject is introduced to the subject pool. Similar improvements are possible that may hold potential clinical application.
Data fusion in human subjects is a challenging problem. Unlike fusion in other areas of imaging, in human brain imaging within subject and between subject inconsistencies are unknown or complicated to compute. It is difficult to conceptualize a correct way to fuse data without making many assumptions. In this paper we introduced straightforward approaches with minimal assumptions to correlate structure with function incorporating all brain voxels and demonstrate methods by which this fusion can reveal differences between patients with schizophrenia and healthy controls. To our knowledge a fusion analysis of this nature is the first of its kind.
We introduce methods to fuse structural and functional imaging data by including all possible combinations of correlations from the whole brain. We show how these methods can be used to extract features that differentiate a group of patients with schizophrenia from a group of healthy controls. Our findings indicate that the correlation between structural data, from gray matter concentration, and functional data, from a sensory motor task activation map, are weaker in patients than in controls. Structural regions in the cerebellum had higher positive correlations with functional data in controls than in patients. Structural regions in the frontal regions had higher negative correlations with functional data in controls than in patients. The results reported in this paper reveal interesting findings that are not possible to derive from conventional one dimensional/unimodal approaches.
The data collection was funded by the Department of Energy, grant DE-FG02-99ER62764. The authors would like to thank the University of Iowa Hospital, Massachusetts General Hospital, the University of Minnesota and the Mind Research Network staff for their efforts during data collection, preprocessing and analysis. This work was supported by the National Institutes of Health under grants 1 R01 EB 006841 and 1 R01 EB 005846.