|Home | About | Journals | Submit | Contact Us | Français|
Developments in multi-channel radio-frequency (RF) coil array technology have enabled functional magnetic resonance imaging (fMRI) with higher degrees of spatial and temporal resolution. While modest improvement in temporal acceleration has been achieved by increasing the number of RF coils, in parallel data acquisition techniques, the maximum attainable acceleration is intrinsically limited only by the amount of independent spatial information in the combined array channels. Since the geometric configuration of a large-n MRI head coil array is similar to that used in EEG electrode or MEG SQUID sensor arrays, the source localization algorithms used in MEG or EEG source imaging can be extended to also process MRI coil array data, resulting in greatly improved temporal resolution by minimizing k-space traversal during signal acquisition. Using a novel approach, we acquire multi-channel MRI head coil array data and then apply inverse reconstruction methods to obtain volumetric fMRI estimates of blood oxygenation level dependent (BOLD) contrast at unprecedented whole-brain acquisition rates of 100 ms per sample. We call this combination of techniques magnetic resonance Inverse Imaging (InI), a method that provides estimates of dynamic spatially-resolved signal change that can be used to construct statistical maps of task-related brain activity. We demonstrate the sensitivity and inter-subject reliability of volumetric InI using an event-related design to probe the hemodynamic signal modulations in primary visual cortex. Robust results from both single subject and group analyses demonstrate the sensitivity and feasibility of using volumetric InI in high temporal resolution investigations of human brain function.
Hemodynamically based fMRI (Belliveau et al., 1991; Belliveau et al., 1990; Kwong et al., 1992; Ogawa et al., 1990) is typically limited to a temporal sampling period of two to four seconds if whole brain coverage is desired. Most fMRI data acquisition methods employ an EPI technique that utilizes many phase encoding steps and multiple read-out gradients. Consequently, this reliance on gradient encoding results in long image acquisition times and relatively loud acoustic noise related to the requisite rapid gradient switching. Here, we demonstrate the use of a novel volumetric imaging method, called Inverse Imaging (InI), which uses minimal phase encoding to achieve an order-of-magnitude improvement in BOLD-contrast temporal resolution. Its minimal dependence on encoding gradients allows extremely short image acquisition times, with an associated trade-off involving somewhat reduced and spatially-varying spatial resolution.
The temporal resolution of MRI is limited by the time required to traverse k-space during signal acquisition. The collection of volumetric MRI data continues until the completion of k-space traversal in multiple 2D k-spaces or in a single 3D k-space. Classical gradient-echo or spin-echo image acquisition methods collect data from one k-space line during each excitation. Thus the total acquisition time for traditional 3D MRI data acquisition is the product of the number of slices and the number of phase encoding steps. In contrast to gradient-echo or spin-echo imaging, both echo-planar imaging (EPI) (Mansfield, 1977) and spiral imaging (Blum et al., 1987) utilize fast gradient switching to achieve 2D k-space traversal in a single RF excitation. With current state-of-the-art EPI or spiral imaging techniques, one 2D single slice image can be collected in approximately 80 ms, allowing whole head coverage with 3 mm isotropic resolution in two to four seconds. Small improvements in temporal resolution can be achieved by optimizing k-space sampling schemes and reconstruction methods: e.g., instead of completing the k-space traversal for every measurement, MRI data acquisition can be accelerated by coordinated alterations of in k-space trajectories and their associated image reconstruction algorithms, as in partial-k space sampling (McGibney et al., 1993). Alternatively, a priori information-based methods can improve the temporal resolution of MR dynamic measurements (Tsao et al., 2001).
Recently, parallel imaging methods have been introduced to reconstruct images using spatial information derived simultaneously from multiple coil array channels The techniques employed include k-space SMASH (Sodickson and Manning, 1997), k-space GRAPPA (Griswold et al., 2002) and image domain SENSE (Pruessmann et al., 1999), all of which share a similar theoretical background (Sodickson and McKenzie, 2001). While parallel MRI can accelerate data acquisition rates by reducing total k-space traversal at the cost of reduced signal-to-noise ratio (SNR), the resulting net acceleration rate is limited both by the number of array coils and the specific phase-encoding scheme employed.
Prior information can be incorporated by combining EPI with parallel MR imaging (Golay et al., 2000; Preibisch et al., 2003; Schmidt et al., 2005; Weiger et al., 2002), resulting in fMRI detection sensitivity improvements with sensitivity encoded parallel MRI techniques (Lin et al., 2005b). Prior-informed parallel MRI has been explored using a fixed regularization parameter with empirical singular value decomposition truncation (King, 2001; Sodickson, 2000). Incorporation of prior information can suppress noise amplification in parallel MRI reconstruction (Lin et al., 2002; Lin et al., 2005b; Lin et al., 2004; Tsao et al., 2002) and traditional parallel MRI has been used to solve under-determined ill-posed problems (Katscher and Manke, 2002; Tsao et al., 2003). However, either only minor acceleration has been achieved (4-fold acceleration using a 2-channel array in cardiac imaging) (Katscher and Manke, 2002), or the reconstruction process has depended on incorporation of low-resolution prior image information (Tsao et al., 2003; Tsao et al., 2005).
More extreme accelerations in MRI acquisition rates have been achieved by reconstructing each image from a single echo. For example, single-echo-acquisition (SEA) was achieved using a dedicated 64-channel linear planar array that eliminated phase encoding, instead using the spatial information obtained from an array of long and parallel coils. This planar pair element design proved to be crucial for achieving well-localized field sensitivity patterns (McDougall and Wright, 2005). In other work, Hennig developed the one-voxel-one-coil (OVOC) MR-encephalography technique, obtaining a reconstructed image by computing the product of a full FOV reference scan and the accelerated acquisition scan where traditional phase and frequency encoding can be selectively omitted. This approach uses simultaneous multi-channel acquisition with multiple small receiver coils sampled such that the signal received by each coil is read out separately. The effective voxel size observed by each receiver channel is determined by the sensitive volume of the corresponding coil element and the source spatial distribution is estimated by constrained reconstruction using images from each separate coil as references (Hennig et al., 2007). A similar reconstruction algorithm termed HYPR was also proposed in the context of MR angiography (Mistretta et al., 2006). Nevertheless, none of these approaches explicitly formulate the relationship between the spatial information contained in the different channels of a RF coil array with full gradient encoding or with minimal gradient encoding. Nor do they provide algorithms to estimate the significance of task-related signal changes that would allow dynamic statistical inferences to be made from a highly temporally resolved data set.
Mathematically, the maximum acceleration possible with parallel MRI acquisition is limited by the available independent spatial information encoded by the coil array elements. This limit manifests itself as a problem in solving an over-determined linear system. Increasing the number of channels can thus increase MRI sampling rates. To this end, dense head coil arrays consisting of 16 (Bodurka et al., 2004; de Zwart et al., 2002; de Zwart et al., 2004), 23, 32 and 90 elements (Wiggins et al., 2005a; Wiggins et al., 2005b) have been constructed in support of a range of parallel acquisition applications. In addition, a dedicated 64-channel linear planar array has been developed to achieve 64-fold acceleration (McDougall and Wright, 2005). Notably, the geometric configuration of our 32-channel head array is remarkably similar to that used for electrode and super-conducting quantum interference (SQUID) sensor arrays in modern EEG and MEG systems (Hamalainen et al., 1993). While the MEG sensors detect magnetic fields generated by neural currents (Hamalainen et al., 1993), MRI detects oscillating electromagnetic fields from magnetization precession (Haacke, 1999). In addition, while MEG derives all of its spatial information from the geometry of the detectors, current accelerated MRI methods still rely heavily on gradient encoding.
We have generalized parallel MRI reconstruction techniques to exceed the limitations encountered when utilizing an under-determined linear system by introducing single-shot volumetric MR Inverse Imaging (InI), an approach that employs an over-determined linear system in order to achieve dramatically reduced acquisition times. We demonstrate the use of single-shot volumetric InI in supporting dynamic spatially-resolved statistical inference in a functional neuroimaging experiment. Inspired by MEG and EEG source localization techniques, we use a generalization of prior-informed parallel MRI (Lin et al., 2005b; Lin et al., 2004) and an adaptation of MEG reconstruction methods to MRI, to reduce the whole-brain sampling time by minimizing the k-space traversal time. Rather than relying on gradient encoding, InI derives spatial information by solving the inverse problem utilizing information from all array channels. Thus, given the constraint imposed by the need to use echo times (TE) that are optimal for BOLD-contrast, InI can complete k-space traversal and acquire sufficient data for whole-brain image reconstruction in under 100 ms. Although we have previously shown the feasibility of a 2D InI implementation (Lin et al., 2006), we now demonstrate its application to functional imaging studies employing 3D whole-brain coverage and event-related designs (Rosen et al., 1998). Event-related fMRI is a widely utilized neuroimaging method to study not only spatial but also temporal activity of hemodynamic changes secondary to the neuronal events. Compared to the classical block-design fMRI, the timing information available in the event-related fMRI allows for the study of both transient and steady states of cerebrovascular responses. This experimental technique mitigates the difficulty of potential bias originated from the contexts or the history of previous stimuli events. Event-related fMRI also enables the analysis of data using post-hoc categorization (Wagner et al., 1998). In some experimental designs, such as “odd-ball” experiments, can only be implemented using event-related fMRI rather than block-design (Friston, 2007). All reasons described above encouraged us to study the feasibility of 3D InI acquisitions and reconstructions using event-related fMRI design.
The principal novelty of our method is its combination of dense coil arrays with a linear estimation approach, allowing the transition from a largely gradient encoded to a largely detector encoded image, thereby achieving an order-of-magnitude speedup in the frame rate of dynamic whole-brain MRI. In the following sections, we present the data acquisition strategy and mathematical algorithms underlying InI, which allow extension of the technique to include event-related, 3D functional imaging designs with whole-brain coverage. We next demonstrate the technique’s capabilities in measuring the spatiotemporal properties of the hemodynamic response to brief visual stimuli using a 100 ms temporal sampling rate on a 3T scanner with a 32-channel coil array. In comparison with conventional EPI, 3D InI exhibits comparable sensitivity and adequate spatial resolving power in detecting visual task-related activity when performance is examined at both the single subject and group levels.
Five healthy participants were recruited for the study. Informed consent approved by the Institutional Review Board was obtained from each participant.
Our participants were asked to maintain fixation at the center of a tangent screen while viewing a high-contrast visual checkerboard reversing at 8 Hz. The checkerboard subtended 20° of visual angle and was generated from 24 evenly distributed wedges (15° each) and eight concentric rings of equal width. The stimuli were generated using the Psychtoolbox (Brainard, 1997; Pelli, 1997) and presented using Matlab (Mathworks, Natick, MA). The checkerboard stimuli were presented for 500 ms duration and the onset of each checkerboard presentation epoch was randomized with a uniform distribution of inter-stimulus intervals varying from 3 to 16 s. Thirty-two visual stimulation epochs were presented during each of four 240 s runs, resulting in a total of 128 visual stimulation epochs per participant.
MRI data were collected on a 3T MRI scanner (Tim Trio, Siemens Medical Solutions, Erlangen, Germany). We used a custom-built 32-channel head array receive coil (Wiggins et al., 2006) and a body transmit coil. The array consisted of 32 circular surface coils tessellated evenly to cover the whole brain. Functional imaging included standard BOLD-contrast imaging using conventional EPI techniques as well as volumetric InI.
Using 3D InI, each acquired volume was obtained at a particular time instant by the combination of EPI frequency encoding along the inferior-superior direction and phase encoding along the anterior-posterior direction. Figure 1 shows a schematic diagram of the spatial encoding procedure used to combine EPI and InI. The spatial information in the left-right direction was recovered during the InI reconstruction computation.
InI reconstruction requires collection of a reference scan that provides information about the entire 3D volume. With this reference scan data, also called the forward operator, accelerated acquisition is enabled by replacing time consuming spatial encoding, dependent upon gradient switching, with an alternative approach utilizing an image reconstruction algorithm.
The reference scan was implemented using a single-slice echo-planar imaging (EPI) readout approach. Specifically, we excited one thick sagittal slab covering the entire brain (FOV 256 mm × 256 mm × 256 mm; 64 × 64 × 64 image matrix), setting the flip angle to the Ernst angle of 30°. Partition phase encoding was used to obtain the spatial information along the left-right axis (inter-aural line). The EPI readout had frequency and phase encoding along the superior-inferior and anterior-posterior axes respectively. We used 100 ms TR, 30 ms TE, 2604 Hz bandwidth and a 10 s total acquisition time for the reference scan.
InI functional data acquisition used the same volume prescription, TR, TE, flip angle, and bandwidth as the reference scan. The principal difference was that the partition phase encoding was removed. The full volume was excited like in the reference scan, and the spins were spatially encoded by a single-slice EPI acquisition. This resulted in a projection image along the left-right direction. The InI reconstruction algorithm, described in the next section, was then used to estimate the spatial information along the x–y axis. In each run, we collected 2,400 measurements after collecting 32 measurements in order to reach the longitudinal magnetization steady state. A total of four runs of data were acquired from each participant.
To validate the InI functional results, conventional EPI data, using conventional BOLD-contrast detection methods, was also collected using identical stimulus and presentation paradigm timing. EPI functional data acquisition used TE = 30 ms, TR = 2 s, flip angle = 90°, FOV = 220 mm, 24 slices, slice thickness = 4 mm with 20% gap).
In addition to the EPI and InI functional data, anatomical MRI data for each participant were obtained in the same session using a high-resolution T1-weighted 3D sequence (MPRAGE, TR/TE/flip = 2530 ms/3.49 ms/7°, partition thickness = 1.33 mm, matrix = 256 × 256, 128 partitions, FOV = 21 cm × 21 cm). Using these data, the location of the gray-white matter boundary was estimated with an automatic segmentation algorithm to yield a triangulated mesh model with approximate 340,000 vertices (Dale et al., 1999; Fischl et al., 2001; Fischl et al., 1999). The spatial registration between the EPI data or InI reconstructed data and the anatomical data was done by using SPM5 (http://www.fil.ion.ucl.ac.uk/spm/), estimating a 12-parameter rigid body transformation between the EPI mean image or volumetric InI reference scan and MPRAGE anatomical study. Note that the actually accelerated InI acquisitions were not spatially registered using this tool. The registration matrix was subsequently applied to each time instant of the EPI or InI hemodynamic estimates, to transform the neural activity estimates for each functional run to an inflated cortical surface space (Dale et al., 1999; Fischl et al., 1999). The transformed results were also spatially smoothed with a 3D Gaussian kernel with 10 mm full-width-half-maximum (FWHM).
Both reference scans and InI scans use EPI read-out, which demonstrate classical Nyquist ghost artifact in the phase-encoding direction (anterior-posterior direction in our implementation). To correct this artifact and to make consistent measurements between reference and InI scans, we performed the same phase correction algorithm in each reference (with partition encoding) and InI (without partition encoding) image at all channels of the coil array respectively. Specifically, 64 interleave EPI read-outs were acquired in the frequency encoding direction (superior-inferior direction) without phase encoding blips to estimate the k-space shift between even and odd echoes along the frequency encoding direction. Such an echo-shifting was then corrected by appending a compensation phase to the even echoes to avoid N/2 Nyquist ghost. The InI acquisition and reference scans were processed using 2D and 3D fast Fourier transformations from the k-space domain to the image domain, respectively. The reference scan in each channel of the coil array was synthetically averaged across partitions to simulate the InI acquisitions by making projection images along the partition encoding direction in each channel of the coil array, such that the simulated InI acquisition at location and channel i , is calculated as
. Here represents the spatial location indices across different partition phase encoding steps with the same frequency and phase encoding numbers, indicated by the spatial index , and Ai () represents the reference scan image from location and channel i . These simulated data were compared with the InI acquisition at each time instant to separately investigate the phase difference between the simulated InI projection image from the reference scan and the actually acquired InI projection image at each time instant. The global phase difference, θi(t) , for channel i at time instant t, is given by
where di (,t) represents the signal from the InI acquisition with spatial location , time t and channel i of the coil array. Phase information is important in the subsequent InI reconstruction since phase may change dramatically over continuous scans as a result of scanner instabilities such as gradient coil heating. We corrected this discrepancy in phase between the actual accelerated InI acquisition and the simulated InI acquisition from the reference scan by subtracting θi(t) for different channel and time instant separately.
After pre-processing of the whole time series in the InI projection acquisitions, now we tried to estimate the hemodynamic response function (HRF) in each projection images across all channels of the coil array. This effort reduces the size of data in time-domain dramatically without affecting the subsequent image-domain image reconstruction. From a list of the stimulus onset times, we constructed a vector with entries containing a one indicate the occurrence of stimuli during each 500 ms stimulation period, and all other entries contain zeros. A contrast matrix D was constructed from the convolution between the vector and the hemodynamic response function H ,
, where denotes the linear convolution.
We used a finite-impulse response (FIR) basis for the HRF and thus H is a discrete delta function with different delay indicating the HRF at a specific time instant. Specifically, H modeled a HRF of 30 second duration with six second pre-stimulus baseline of 100 ms TR. Thus H is an identity matrix of dimension 300 (i.e., nh = 300). In addition to the phase drift θi(t), we appended the design matrix D with two confound vectors within each run to account for linear drift and constant confounds. Additionally, we also included a time series from the mean of the magnitude projection image in each channel as a confound vector attempting to explain the potential physiological noise (respiration and cardic pulse). The estimation of GLM coefficients for channel i and time instant τ of the InI acquisition was computed by least square fitting
where the superscript H indicates the transpose and complex conjugate. HRF estimation, hi(,τ), was then extracted from entries of βi(,τ) vector corresponding to FIR basis function while all other contribution corresponding to confounding vectors was discarded in the subsequent analysis.
The minimum-norm estimate (MNE) provided the reconstructed image x(,τ) of 3D volumetric index and time index τ , which can be expressed as
where C is the noise covariance matrix of the array, λ() is the regularization parameter, and WMNE() is the inverse operator for the InI location index . A() and h(,τ) indicate the vertically concatenated Ai() and hi(,τ) spatial index and time index τ across all channels in the coil array. The regularization parameter was calculated from a pre-defined signal-to-noise-ratio (SNR) as
Here Tr(•) indicates the trace of the matrix. This choice of regularization parameter followed our original 2D InI reconstruction algorithm and was inspired by the analogy between InI reconstruction (Lin et al., 2006) and distributed source modeling in MEG (Lin et al., 2005a). We used SNR = 5 in this study. Such a choice of SNR was derived from our previous InI study (Lin et al., 2006). Sensitivity of MNE reconstruction in InI was investigated in our previous study, which showed that the regularization parameter affects InI reconstructions moderately if it changes within ±10 fold. Our choice of regularization parameter should be within that range.
To allow statistical inference from the results of the InI time series reconstruction, we estimated the noise in the baseline by applying the MNE inverse operator to the baseline InI data. Then, dynamic statistical parametric maps were derived as the time-point by time-point ratio between InI reconstruction and this baseline noise estimates, given by
where diag(•) is the operator to construct a diagonal matrix from the input argument vector. Here x(,τ) represents estimated signal vector and diag(WMNE()CWMINE()H) denotes estimated noise vector. Both vectors have the dimension of the number of sources to be resolved. The division denotes the element-wise division. Dynamic statistical parametric maps (dSPMs) t(,τ) should be t distributed under the null hypothesis of no hemodynamic response (i.e., x(,τ) = 0) (Dale et al., 2000). When the number of time samples to calculate the noise covariance matrix C is quite large, the t distribution approaches the unit normal distribution (i.e., a z-score).
The statistical modeling of the conventional EPI BOLD-contrast data was implemented using a general linear model (GLM) with a design matrix based on the same finite impulse response (FIR) model (c.f. Eq. 3 and Eq. 4). Using least square fitting, we estimated the strength of the FIR HDR at each voxel at each voxel and each time instant, and the residuals in GLM led to an estimate of baseline variability. Taken together, the t-statistics of the evoked hemodynamic response were calculated for each image voxel at each time instant.
Here we want to emphasize that InI image reconstruction had separate processing for HRF estimation in the projection images in all channels of the coil array (time domain) and restoration of 3D spatial information using MNE (spatial domain). However, EPI data had only time domain HRF estimation. In the time domain estimation of HRF, both InI and EPI data used the same finite-impulse-response (FIR) basis in the GLM estimation of HDR and thus the results are comparable. The use of FIR basis function has no prior assumption on the shape of HDR. This is different from the traditional use of parameterized models to estimate canonical HDR using, for example, Gamma functions. Compared to traditional EPI reconstruction, the only additional processing in InI reconstruction is the MNE in the spatial domain to convert the 2D projected images into 3D volumetric images.
To compare the similarity of InI and EPI reconstructions, InI dSPMs of t-statistic maps were averaged between 3–7 s after stimulus onset. The center of mass of this temporally averaged t-statistic map from the InI TINI() and EPI t-statistics maps TINI() were calculated separately:
Here ‖•‖ denotes the Euclidean norm.
We performed numerical simulations to evaluate the spatial resolution and localization accuracy of our InI reconstructions. The reference data for the forward operator A() and noise covariance matrix C . The simulation procedure started from creating a source vector x(), with asset to unit activity and other locations set to zero. We then estimated the idealized measurements from all coil array channels by computing the product of the forward operator A() and x().
We created 100 realizations of synthetic noise with spatial coloring according to the noise covariance matrix:
where ε() is the noise vector with complex values following a Gaussian distribution of zero mean and unit variance. Uc and Sc are the singular vectors and singular values of the noise covariance matrix. At a specified SNR, the noise n() was scaled and subsequently added to s() to generate the synthetic measurements d():
Then we followed Eq. 5 to obtain the data covariance matrix, the inverse operator, and the noise normalized inverse operator:
. The InI reconstruction obtained with this procedure is equivalent to the point spread function of the simulated source x() . Both were scaled to a maximum of 1.
Similar to procedures used in MEG/EEG source analysis (Dale et al., 2000; Liu et al., 1998; Liu et al., 2002), we estimated the average point spread function (aPSF) at each location to quantify the spatial spread of the reconstruction:
where |di()| indicates the distance between source location i and source location . represents vector entries in the InI reconstruction exceeding 0.5 and l is the number of voxels to be spatially resolved by the InI reconstructions. This procedure allows estimation of the full-width-half-maximum (FWHM) of the point spread function. A 3D map of the spatial distribution of the average point spread function, aPSF, for either MNE or MNE dSPM estimates can be obtained by repeating the calculation across the whole source space, the 256 × 256 × 256 mm FOV and 4 mm isotropic spatial resolution.
Since InI is an intrinsically ill-posed inverse problem, the reconstructed image may not be the original spatial distribution of spins contributing to the actual measurements. Thus the analysis of localization accuracy on the discrepancy between the reconstructed sources and the original sources is desired. Quantification of localization accuracy was done by calculating the shift between the center of mass of InI reconstruction and the simulated source:
A 3D SHIFT metric map for MNE and MNE-dSPM was generated in each simulation. Since the inverse operators depend on the SNR (and the measurement data), the SNRs were parametrically varied from 0.1 to 100. Figure 2 illustrates the procedure of spatial resolution analysis using simulations. The implementation of the image reconstruction and statistical analysis procedures were done with Matlab.
Figure 3 shows the spatial distribution of aPSF at different SNRs (SNR = 0.5, 1, 5, and 10). We observed global reductions in the aPSF for MNE reconstructions at higher SNRs. In particular, deep brain regions show a larger aPSF metric. This matches the physical intuition that at center of the head coil the B1 fields from all channels are less spatially disparate as well as the SNR is lower than that at periphery and thus the spatial resolution becomes worse. Quantitatively, the MNE-dSPM inverse still has an average aPSF of 8.64 mm when SNR is higher than 1. At extremely high SNR (SNR >50), MNE can provide an excellent spatial resolution (average aPSF of 0.15 mm). To translate the aPSF into a measure of spatial resolution, the point spread function reported here should be spatially convolved with the nominal spatial resolution of the fully gradient encoded scan. Thus, for example, the average spatial resolution at SNR=5 is approximately 8.7 mm with a standard deviation of 6.5 mm (see Table 1).
The SHIFT metrics derived from MNE were shown in Figure 4. MNE is characterized by sporadic high SHIFT metrics at low SNRs (SNR =0.5 and 1). Larger errors occurred at source locations in deep brain areas. On average, the localization accuracy is higher than 5 mm when SNR is higher than 1. Details of the aPSF and SHIFT metrics are listed in Table 1.
Raw traces of InI acquisitions from three channels close to the occipital lobe, parietal lobe, and frontal lobe were shown in Figure 5. Clear cardiac and respiratory fluctuations were observed in the phase and magnitude plots of the time series from each RF coil channel. Figure 5 also shows that acquired InI images indeed had significant phase drifting at respiratory fluctuation. This supports the requirement of phase correction in Eq. 2. To illustrate the quality and consistency of the InI reconstruction, we show the difference magnitude images between simulated InI acquisitions from reference scan and the those images of the mean of actually accelerated InI acquisitions in Figure 6. The discrepancy between the simulated and actual InI acquisitions was quantified as the percentage error with respect to the simulated InI acquisitions in each channel of the coil array separately. We observed that all 32 channels show a discrepancy of less than 10%. Note that the number of channel is arbitrary and it does not imply spatial location of each RF coil. This implies that the reference scan suffices an accurate forward operator for the subsequent inverse operator derivation and volumetric image reconstruction.
To illustrate the spatial resolution of InI reconstruction, we performed InI reconstruction using the average InI acquisitions across time. The results are shown in Figure 7, which includes magnitude images of InI reconstructions and sum-of-squares reference images across 32-channels of the coil array at central 12 contiguous axial slices. Compared to the reference images, InI reconstructions preserved some features, including boundary, cerebral hemisphere boundary, and the contours of the whole brain. However, localized image features and contrasts are different. In the context of fMRI, we are interested in the time-domain contrast-to-noise ratio, rather than the spatial-domain contrast-to-noise ratio. Also, using univariate general linear model, the analysis treats different image voxel separately. Thus we tolerated such spatial heterogeneity and proceeded with time-series analysis.
Figure 8 shows two series (2.0 to 2.9 s and 6.0 – 6.9 s post stimulus) of functional activity estimated by InI dSPM overlaid on the sum-of-squares images from the reference scans at eight axial slices from a representative participant. The InI dSPM t-statistic maps were rendered on an inflated brain surface where light gray indicates gyri and dark gray indicates sulci. We used a critical threshold of t = 2 (uncorrected p-value = 0.028). The maps show a progressive activation in response to the reversing checkerboard around the calcarine sulcus. Starting from the most posterior part of the occipital lobe, the activity spreads anteriorly, superiorly, and inferiorly. Peak activity was observed between three and five s after stimulus onset. This response started to decrease approximately six seconds after stimulus onset
The area showing a positive visual response in the first three seconds after stimulus onset was used to define a region of interest (ROI) in the primary visual cortex (V1). The average values and standard deviations of this ROI time course are shown in Figure 9. Without utilizing any specific model of the hemodynamic response, InI measures revealed a sharp BOLD-contrast signal peak at four s after stimulus onset. We also observed a post-stimulus undershoot between 10 and 24 s after stimulus onset. The pre-stimulus interval shows fluctuations around the baseline.
To compare InI with conventional EPI reconstruction, we compared data from the same participants studied with both techniques. Figure 10 shows the t-statistic maps from both InI and EPI acquisitions using an inflated cortical surface model. We averaged the InI dSPM t-statistics between three and five second after stimulus onset in order to match the two-second TR of the EPI acquisition from the same V1 ROI. To compensate for the differing sensitivity of the two methods, critical thresholds were chosen as t < 2 (uncorrected p-value <0.028) and t < 4 (uncorrected p-value <10−4) for InI and EPI respectively in order to show similar size and pattern of cortical activity.
Differing detection sensitivity between InI and EPI were quantified by examining the t-value maxima, the t-value baseline standard deviation, and the z-score transform of the t-values in primary visual cortex of each participant (see Table 2). Generally, EPI was associated with higher peak t-values compared to InI. However, the baseline variability of the EPI t-values was also higher than InI. The overall InI z-score transform of the t -values is thus higher than EPI. We noticed that the standard deviation of baseline t-statistics does not equal to 1 in V1 ROI. The reason for such a result is because we calculated the averaged t-statistics within V1 ROI first and then we calculated the standard deviation. The averaging within the V1 ROI certainly improves the noise fluctuation because of averaging across voxels within the ROI. This explained why the reported standard deviation of the t-statistics in the baseline period is less than 1.
The distances η between the centers of mass of the EPI and InI reconstructions from five participants are listed in Table 3. On average, the spatial localization of InI and EPI peak activities differs by 5.2 mm. This localization discrepancy was quite stable among all five participants. The maximum distance is 6.8 mm and minimal distance is 3.4 mm, both of which are within 2 mm difference compared to the average.
Figure 11 shows single frames of InI dSPM t-values averaged over five participants at 100 ms temporal resolution. The shape of this group average time series is very similar in character to the individual participant time series shown in Figure 9, demonstrating the stable performance of InI acquisition and reconstruction techniques. Notably, the negative t-values observed in the single participant, shown in Figure 9, disappeared in this group average, implying that the signal decrease, as evidenced by the negative t-values, may be an individual or temporally variable effect. The individual frames of this group average show progressive increasing activity starting at 2.7 s after the stimulus onset (critical threshold t < 2; uncorrected p-value <0.028). The signal returned to baseline approximately 6.0 s after stimulus onset.
The time course of the InI dSPM t-values from the group average are shown in Figure 12. The shape of the hemodynamic response in the group average is similar to that observed in individual participants (see Figure 3). We observed a reduced variability along the time course, potentially resulting from averaging multiple participants. The peak activity was found at 4.5 s after stimulus onset and the post-stimulus undershoot was observed between 10 s and 24 s after stimulus onset in the group average.
We have shown that single-shot volumetric InI methods can achieve an order-of-magnitude acceleration in hemodynamic response sampling by combining dense head coil array parallel data acquisition with distributed source modeling. Applying volumetric InI to study the visual system using an event-related design, we found that the method is both sensitive and reliable, as demonstrated in the individual participant and group average results. Localization accuracy was examined by comparing the InI and conventional EPI task-related activity in primary visual cortex, with the comparison revealing a good spatial match (approximately 5 mm) between the two methods. Our results in the visual cortex activation show faster BOLD response compared to conventional EPI. This is comparable to results from literature (Janz et al., 1997; Pfeuffer et al., 2002). In comparison to previous studies, InI methods exhibit two principal advantages. First, they allow collection of whole-brain volumetric data using single-shot EPI acquisition, while previously it was only possible to obtain one single 2D image at high temporal resolution (20 ms). Volumetric acquisition is more attractive to neuroscientists primarily interested in investigating whole-brain spatiotemporal activity patterns. By sampling the entire brain, volumetric InI can avoid the need for tedious manual slice prescription based on the prior anatomical scans required to identify target brain areas. Second, our approach combines volumetric InI and event-related fMRI design (Rosen et al., 1998), a commonly employed experimental approach in functional neuroimaging studies. We found that volumetric InI acquisition and reconstruction had sufficient sensitivity to reliably detect task-related activity in primary visual areas with high temporal resolution. We also found that this study we had reduced false positive activation compared to our first attempt of InI (Lin et al., 2006) (see Figure 8). This is because the current implementation has better matching between the reference scan and the InI scan: both scan used the same 3D RF excitation, frequency encoding, and phase encoding. The only difference is that the reference scan had partition encoding steps, but InI scans did not. The improved matching between the reference scans and the InI acquisitions, as quantified in Figure 6, helps to generate a better reconstruction with higher sensitivity and reduced false positive activation.
The difference of SNR and sensitivity of EPI and InI have been analytically and empirically reported in our results. Both methods used BOLD contrast. However, due to different image reconstruction algorithm and sampling rate, EPI and InI have different sensitivity and degree of freedom. Compared to EPI, InI has an additional step of solving an inverse problem. Both EPI and InI time series modeling and statistical inference used the same general linear model. Using an MNE inverse operator in InI reconstruction resulted in a “decreased” sensitivity (as reported in the reduced peak t-statistics) as well as “decreased” detection variability (as reported in the reduced baseline t-statistics variation). This is due to the spatial smoothing intrinsic to the minimum L-2 norm inverse operator. Such a smoothing effect was analyzed quantitatively in the average point spread function simulations Note that InI and EPI also have different degrees of freedom, since they are of different sampling rates. We are aware of the potential temporal correlation inside the acquisitions. However, a further detail analysis on such correlation is beyond the scope of this study and we will pursue it in the near future.
One limitation in BOLD-contrast fMRI detection sensitivity is physiological noise originating from cardiac and respiratory fluctuations (Kruger and Glover, 2001). It has been shown that physiological noise sources are the dominant limiting factor in high-field fMRI (Giove et al.). Traditional EPI takes approximately two to four seconds to acquire a full 3D brain volume and therefore lacks sufficient temporal resolution to resolve the physiological noise which then becomes aliased into the recorded time series. It will be possible to utilize the high temporal resolution of InI to improve the sensitivity of BOLD-contrast fMRI by acquiring single-shot multiple-echo datasets at 10 Hz sampling rates. Due to this rapid sampling rate, volumetric InI data can satisfy the Nyquist sampling criterion and thereby allowing more effective digital filtering strategies to isolate and reduce cardiac and respiratory sources of noise. These reductions in temporal noise will result in improved detection of the relevant effects of interest in the experimental design.
InI solves an ill-posed inverse problem in image reconstruction. The limited spatial resolution from all channels of the coil array may not be able to provide a stable and unique solution for image reconstruction. This deficiency in independent coil information leads to limits in spatial resolution. Among all inverse solution alternatives, we want to bring a caution on the nominal spatial resolution and the actual image resolution. For example, using a classical equivalent current dipole (ECD), the spatial resolution can be regarded as infinitely high. However, this does not adequately reflect the physical resolution set by the measurement. The reconstructed InI images using MNE algorithm preserved some low-resolution image figures (Figure 7). This matches to the well-known spatial smoothing effect in a minimum- norm solution. Previously, we estimated that 2D InI usingin MNE reconstruction has a spatial resolution ranging between 5 mm and 16 mm, depending on the measurement SNR and the coil geometry (Lin et al., 2006). From our results, we can see that volumetric InI provides reasonable spatial accuracy relative to conventional EPI acquisition. However, a more systematic and extensive study of all factors modulating the spatial resolution, including SNR, coil geometry, field strength, and regularization parameters, will be required in order to quantify the factors influencing spatial variations in resolving power.
The spatial resolution is spatially-varying because the amount of independent spatial information from all channels in a coil array is varies among the image voxels within the field-of-view. Notably, in peripheral cortical areas the spatial resolution is higher and at the center of the brain the spatial resolution is lower (Lin et al., 2006). In volumetric InI, anisotropic spatial resolution effects only appear in the InI dimension (L-R direction in this study), while the other two spatial dimensions (A-P and S-I directions) still retain isotropic spatial resolution because gradient encoding is used. To improve the spatial resolution, there are two alternatives. First, increasing the number of channels in a coil array can provide more independent spatial information. However, the benefit of increasing channels will reach a plateau as the consequence of electromagnetic theoretical limitations (Ohliger et al., 2003; Wiesinger et al., 2004). In addition, at high field more independent spatial information can be obtained from the same geometry of a coil array as the consequence of a shorter wavelength. This implies that higher spatial resolution can be obtained at 7T, using the existing 32-channel coil array geometry and MNE reconstruction algorithm. Second, we can ameliorate the spatial blurring using other inverse reconstruction kernels, as discussed in the following section.
The volumetric InI technique solves an ill-posed inverse problem to obtain reconstructed images, using the MNE solution to estimate the spatial distribution of task-related activity. In MEG and EEG distributed source modeling, spatial filters using linear constraint minimal variance (LCMV) beamformers have been also investigated extensively (Hillebrand and Barnes, 2005; Robinson, 2004; Sekihara et al., 2002; Van Veen et al., 1997). In contrast to MNE, the LCMV approach minimizes the point spread function of the inverse operator, resulting in LCMV results that are more focal compared to MNE. However, one disadvantage of LCMV is its inefficient detection and separation of coherent sources (Van Veen et al., 1997). Interestingly, it has been shown that MNE and LCMV are mathematically related to each other; the difference residing in the fact that MNE is model-driven while LCMV is data-driven during the estimation of the data covariance matrix (Mosher et al., 2003).
Under the interpretation of Bayesian estimation theory, both LCMV and MNE use L-2 norm prior models to solve an inverse problem. Thus it may be possible to replace the L-2 norm model by an L-1 norm model to achieve a narrower point spread function in InI image reconstruction. The reason that the L-1 norm model has a higher spatial resolution (a narrower point spread function) is that the probabilistic distribution of L-1 norm follows a bi-exponential distribution and the L-2 norm model follows a normal distribution. Moreover, the bi-exponential distribution has a more concentrated probability around the mean compared to the normal distribution (Uutela et al., 1999). L-1 norm models have been previously studied in the context of MEG and EEG source localization under the term of “minimum-current estimate” MCE (Lin et al., 2005a; Matsuura and Okabe, 1995; Uutela et al., 1999). In future work, we will investigate the comparative spatial localization accuracy of the LCMV inverse and MCE approaches using InI data.
In this study we employed a single-shot volumetric InI acquisition with an EPI readout. Thus, the reconstructed projection image contains the expected EPI artifacts, including intra-voxel signal loss due to spatially inhomogeneous susceptibility distribution and geometrical distortion along the phase encoding direction due to its intrinsically slower bandwidth. Correction of these artifacts has been studied extensively. For example, to mitigate these artifacts, we can use field mapping to investigate the spatial distribution of the off-resonance effects and then use this information to reduce the susceptibility artifacts (Chen et al., 2006; Chen and Wyrwicz, 1999; Zeng and Constable, 2002). It is also possible to use parallel imaging techniques with EPI acquisitions to limit geometric distortion (Weiger et al., 2002), by systematically skipping multiple integer lines in the continuous sampling of different phase-encoding lines and then reconstructing the skipped phase encoding lines using spatial information embedded inside different array channels. Thus the effective bandwidth in the reconstructed image will be wider and will thereby reduce the distortion. Note that here the spatial information from the coil array channels is in an orthogonal direction between EPI phase encoding (anterior-posterior axis) and InI encoding (left-right axis). Thus the implementation of this approach will not reduce the available spatial information in the InI reconstructions, while SNR loss is a price that must be paid for reduced data sample numbers and changes in the parallel MRI reconstruction geometrical factor (g-factor).
The brain is a highly dynamic system. Consequently, detection and localization of static patterns of regional functional specialization will not be sufficient to fully understand the neural mechanisms underlying complex behavior. Nevertheless, current fMRI research mainly relies on methods that are most sensitive to these stationary patterns of task-related activity. Volumetric InI allows functional imaging measurements at high sampling rates over extended brain volumes using minimal gradient encoding. To our knowledge, this is the most rapid whole-brain fMRI achieved to date and this relatively high sampling rate may, in future work, enable the measurement of relative brain activity onset times and thereby provide a better understanding of the dynamic driving relationships among interacting neural subsystems across the whole brain. This goal cannot be achieved with conventional EPI methods due to their limited spatial coverage and temporal resolution.
Because of the relative sluggishness of the fMRI BOLD-contrast time course, it has been argued that improving the temporal resolution of fMRI might not reveal more information about brain dynamics. However, it is becoming increasingly clear that this popularly held view is incorrect, possibly resulting from a paucity of practical methods allowing ultrafast fMRI recordings in large brain volumes. For example, using limited spatial coverage, high temporal resolution fMRI (TR = 100 ms per slice or 400 ms per volume) can detect spatially distinct differences in neural activity onset times that can then be correlated with reaction time in order to localize the cortical areas active in the different phases of sensorimotor integration (Menon et al., 1998). In addition, fMRI experiments using jittered inter-stimulus-interval manipulations have detected millisecond-scale interactions among neural systems (Ogawa et al., 2000). However, no direct measurements of BOLD-contrast signals at high temporal resolution with spatial coverage of the entire brain have been achieved previously. Using volumetric InI, it may be possible to adapt time-resolved imaging techniques to experimental designs probing mutual interactions among a spatially distributed set of participating functionally specialized regions.
Motion correction is a critical process to improve detection power of fMRI. However, Motion correction on InI data is difficult because the acquisitions have only projection images. Such reduction of one dimensional spatial information poses a technical difficulty of utilizing post-processing method to do, for example, rigid body transformation, an intrinsic 3D processing. Thus currently InI may not have the optimal sensitivity due to motion artifacts. A potential method to mitigate this motion artifact issue is using navigator echoes during acquisition to correct motion effects on the flight at a moderate cost of temporal resolution (van der Kouwe et al., 2006). Even though volumetric InI allows dramatic improvement in sampling rates, it is still constrained by the need to use the optimal TE for detection of BOLD-contrast effects (approximately 30 ms at 3T). At higher field strengths, such as 7T, the optimal TE for BOLD-contrast would be 20 ms or less, allowing further acceleration of volumetric InI, possibly to 50 ms whole-brain sampling times. In addition, there are two potential approaches to mitigate the temporal resolution limitations: First, we may use different contrast mechanisms, such as steady-state free precession (SSFP), where TE is usually less than 5 ms (Miller et al., 2003; Miller et al., 2006). However, SSFP contrast is challenging at high field because of its high SAR from RF excitation. Second, we may use an echo-shifting technique to reduce the TR (Golay et al., 2000; Liu et al., 1993), as demonstrated in our previous 2D InI study (Lin et al., 2006). Both alternatives are capable of reducing the sampling time to around 20 ms. The resulting higher temporal resolution may allow study of relative cortical activity timing on an extraordinarily fine time scale, thereby facilitating study of the complex interactions among regionally specialized neural subsystems responsible for the mediation of complex behavior.
This work was supported by National Institutes of Health Grants R01 HD040712, R01 NS037462, P41 RR14075, the Mental Illness and Neuroscience Discovery Institute (MIND).
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.