Spatial normalization, bulk motion correction and coregistration for functional magnetic resonance imaging of the human cervical spinal cord and brainstem

Functional magnetic resonance imaging (fMRI) of the cortex is a powerful tool for neuroscience research, and its use has been extended into the brainstem and spinal cord as well. However, there are significant technical challenges with extrapolating the developments that have been achieved in the cortex to their use in the brainstem and spinal cord. Here, we develop a normalized coordinate system for the cervical spinal cord and brainstem, demonstrating a semiautomated method for spatially normalizing and coregistering fMRI data from these regions. fMRI data from 24 experiments in eight volunteers are normalized and combined to create the first anatomical reference volume, and based on this volume, we define a standardized region-of-interest (ROI) mask, as well as a map of 52 anatomical regions, which can be applied automatically to fMRI results. The normalization is demonstrated to have an accuracy of less than 2 mm in 93% of anatomical test points. The reverse of the normalization procedure is also demonstrated for automatic alignment of the standardized ROI mask and region-label map with fMRI data in its original (unnormalized) format. A reliable method for spatially normalizing fMRI data is essential for analyses of group data and for assessing the effects of spinal cord injury or disease on an individual basis by comparing with results from healthy subjects.


Introduction
Functional magnetic resonance imaging (fMRI) has become one of the most powerful tools for neuroscience research, and yet its use is limited to studies of the cortex, with relatively few exceptions [1][2][3][4]. This means that studies of distributed systems, such as those involved with pain, are limited in their scope and completeness until fMRI can be applied in the entire central nervous system (CNS). It has been demonstrated that important neuronal activity involved with modulation of pain and sensory responses can be identified in the brainstem with fMRI [4][5][6]. Similarly, fMRI of the spinal cord (spinal fMRI) has been demonstrated to show areas of activity in the cervical and lumbar regions with high sensitivity and reliability, in response to thermal, sensory, motor and painful stimuli [3,[7][8][9]. Descending modulation of activity in the cervical spinal cord from brainstem regions, corresponding to emotional factors and focused attention, has also been demonstrated during thermal sensory stimulation (unpublished results). Together, these studies demonstrate the feasibility of fMRI spanning the entire CNS.
Despite the advances in fMRI techniques, there are fundamental problems in extrapolating the technological developments in the cortex to the brainstem and spinal cord [1]. These are caused by the unique challenges posed by fMRI of the brainstem and spinal cord, such as the relatively small cross-sectional dimensions and large rostral-caudal (R/C) extent, the variable curvature of the spine between individuals, the relatively poor local magnetic field homogeneity causing image distortion and low signal intensity in standard T 2 ⁎-weighted fMRI [10] and the problem of motion with each heartbeat [11]. These challenges have required fMRI data acquisition and analysis methods to be adapted and modified for use in the spinal cord [1]. Here, we present the first methods for spatial normalization and nonrigid motion correction of fMRI data spanning the entire brainstem and cervical spinal cord. These methods enable voxel-by-voxel comparisons of results across subjects, as well as automatic labeling of anatomical regions, as is routinely done for essentially all fMRI studies of the brain. This development is an essential step toward obtaining fMRI data that span the entire CNS. It permits group analyses, enables automatic identification of active regions and facilitates comparisons with normative data for assessing injury or disease in the spinal cord or brainstem. fMRI methods that have been developed for use in the spinal cord are a logical choice for fMRI of the brainstem because of the similar challenges that are encountered. Most fMRI studies of the brainstem to date that have succeeded in detecting activity with good reliability have been based on conventional brain fMRI methods [T 2 ⁎-weighted echoplanar imaging (EPI)] [5,6,[12][13][14]. However, spin-echo imaging methods are considerably less sensitive to magnetic field inhomogeneities than gradient-echo methods, and fast spin-echo methods suffer less distortion than spin-echo methods with an echo-planar spatial encoding scheme. As a result, the optimal spinal cord fMRI method is based on a single-shot fast spin-echo [1,15]. The small transverse dimensions and large R/C extent of the spinal cord have been accommodated by imaging thin, contiguous sagittal slices to obtain three-dimensional (3D) functional image data with relatively high resolution [15]. Motion artifacts have been shown to be reduced by means of flowcompensation gradients and spatial saturation pulses to eliminate signal arising from anterior to the spine (i.e., heart, lungs and throat).
The method proposed for the present study is based on an established spinal fMRI method [15,16], with modifications to span the entire cervical spinal cord and brainstem. This choice of imaging method is significant for the proposed normalization procedure because it produces images with high spatial accuracy, unlike methods based on EPI spatial encoding. Here, we define a coordinate system based on the anatomy of the cord and demonstrate a method for normalizing the anatomy to consistent dimensions and posture. The results also enable definition of a standard normalized reference volume and a standardized map of anatomical regions. The normalization procedure can also be applied in reverse, to conform the map of anatomical regions onto the original image data, regardless of spinal cord size or curvature.

Image acquisition
The data used in the present study (for the assessment of the normalization and motion correction methods) are part of a larger study of thermal responses in the brainstem and spinal cord, but the fMRI results themselves are not the topic of this article. Spinal fMRI studies were carried out in a 3-T Siemens Magnetom Trio using a phased-array spine receiver coil with subjects lying supine. Eight healthy volunteers were studied and provided informed consent prior to participating. Initial localizer images were acquired in three planes as a reference for slice positioning for subsequent fMRI studies. The functional image data were acquired with a half-Fourier single-shot fast spin-echo sequence, with the echo time set at 38 ms and with a repetition time of 1 s per slice, to obtain essentially protondensity-weighted images. Signal intensity changes observed upon a change in neuronal activity were the result of signal enhancement by extravascular water protons, as well as a contribution from the BOLD effect, as described previously [17,18]. Sagittal image data were acquired from 14 contiguous sagittal slices, each 2 mm thick, making the effective repetition time equal to 14 s. Each slice was centered roughly on the second cervical (C2) vertebra, spanning from the C7/T1 disc to the top of the thalamus with a 20 cm×10 cm field of view (FOV) and a 192×96 acquisition matrix. The resulting voxel dimensions were therefore 1.04 mm×1.04 mm×2.00 mm. As previously reported, spatial suppression pulses were applied anterior to the spinal canal (to eliminate signal from this region), and flow-compensation gradients were applied in the R/C direction to reduce flow-induced image artifacts.

Motion correction and spatial smoothing
Correction for bulk motion of the spine and spatial smoothing are addressed at the same time as spatial normalization because they are related processes requiring spatial transformations. The methods are applied in sagittal planes to avoid combining data from slices acquired at different times and because the most common bulk motion occurs in the anterior-posterior (A/P) and R/C directions.
Custom-made software written in MatLab (The Mathworks Inc., Natick, MA) was used to first draw a reference line along the anterior edge of the cord in a midsagittal slice. This line was also extended along the brainstem in a straight line, originating from the rostral end of the medulla and running tangent to the anterior edge of the thalamus (Fig. 1) to account for differences in neck flexion and head position. The caudal edge of the pons and the C7/T1 intervertebral disc were marked along this line as position reference points.
The data were then smoothed spatially, but only in the direction parallel to the long axis of the cord (i.e., parallel to the manually defined reference line), in order to avoid introducing partial-volume effects. This was achieved by defining a two-dimensional (2D) kernel to smooth along a line parallel to one axis of the image and then rotating the kernel to be parallel to the closest section of the reference line. The direction of the reference line varied slowly and smoothly across the imaging FOV in every case; hence, smoothing was applied piecewise to each image in 15×15 voxel sections to reduce the computation time.
Bulk motion that may have occurred during the acquisition time series was corrected using functions in the MatLab Image Processing Toolbox. The process was applied on a slice-by-slice basis to align each image in the time series to the first image. First, a set of reference points was defined so as to overlay the entire spine. This included points along two lines parallel to the manually drawn reference line (displaced by 5 pixels anteriorly and 10 pixels posteriorly from the reference line), as well as points on a regular 2D grid (spaced 10 voxels apart in each direction) spanning the entire image. This choice of reference points encompassed the entire image, with a greater concentration of points along the spine, so that the alignment of the spine dominated the image transformation. The "cpcorr" function in MatLab determined the displacement of each reference point to maximize the local correlation of the two images. The polynomial transformation required to optimally align all of the points was then determined with "cp2tform," and the image transformation was applied using the "imtransform" function. The result was nonlinear coalignment of the images that accounts for changes in curvature of the spine over time, as well as translation in the A/P and R/C directions but not in the left-right (L/R) direction.

Spatial normalization
The spatial normalization procedure can be applied as needed to anatomical images, to the results of fMRI data analysis or to each volume of the time-series data to enable group analyses. Applying the normalization after the analysis enables the results to be displayed in both original and normalized formats and eliminates the potential of additional spatial smoothing that may be imposed if the interpolation in the normalization process was applied before the analysis. The normalization procedure consists of aligning the data with a 3D coordinate system with one dimension parallel to the long axis of the cord and the other two running A/P and L/R. The 3D image data were first interpolated to cubic voxels (1 mm in each direction) to make this transformation, and the volume was resliced transverse to the reference line every 1 mm along the R/C direction. A linear transformation was then applied to the volume to place the reference pointsat the caudal edge of the pons (the pontomedullary junction) and at the C7/T1 disc (shown in Fig. 1)at 65 and 205 mm, respectively, from the rostral end of the volume. The span of 140 mm between the two reference points was based on typical dimensions observed in data from a small number of volunteers.
Initially, a reference volume was created by averaging across all time points of a set of spinal fMRI data in which the cord and brainstem were centered in L/R and A/P directions and did not demonstrate any curvature. The data from each individual subject were then aligned with this reference volume by means of rigid-body L/R and A/P shifts at each R/C position, thus fine-tuning the normalization and correcting for any L/R and A/P curvature. A normalized reference volume spanning the brainstem and cervical spinal cord was then constructed by averaging 24 sets of time-series fMRI data from eight healthy volunteers.
Based on the normalized reference volume, a standard region-of-interest (ROI) mask was created. This is a binary mask with values of 1 in pixels overlying the spinal cord or brainstem and 0 otherwise. The mask enables automatic selection of image data from the cord and brainstem and exclusion of surrounding anatomy (CSF, vertebrae, etc.). Based on this ROI mask, a map to label anatomical regions was then constructed by identifying the R/C positions of the thalamus, midbrain, pons, medulla and approximate limits of each cervical spinal cord segment. These R/C regions were then further subdivided into dorsal-ventral and L/R sections, and in the medulla, a midline region was also defined, resulting in a total of 52 anatomical zones. A 3D volume was then constructed with each pixel labeled according to its corresponding zone, resulting in a standardized anatomical region-label mask. The region-label mask and the ROI mask can be used with the normalized data, or the inverse of the transformation used to normalize the image data can be applied to the masks, mapping them onto the original (unnormalized) fMRI data.

Validation of the normalization procedure
The normalization results were assessed quantitatively by determining the error between corresponding anatomical points, a method described by Grachev et al. [19]. Using this approach, a large number of reference points (n=1000) were defined on a regular 2D grid spanning the cord and brainstem, as well as adjacent structures in a midline sagittal slice. Using the MatLab function "cpcorr," the optimal displacement of each reference point (up to 4 pixels along each axis) was determined to maximize the local correlation between the normalized image data and the reference images. The required displacement was taken as the error in the normalized results at each reference point.

Results
The efficacy of the normalization method across the cervical spinal cord and brainstem is shown in Fig. 2, as an averaged image across 24 experiments in the eight volunteers studied. Anatomical details such as the vertebral processes and the well-defined edges of the spinal cord and brainstem structures demonstrate that these features were precisely aligned across the volunteers.
The accuracy assessment, based on anatomical correspondence, showed that the mean error was 0.3 mm and the median was 0 mm. Peak errors were approximately 5 mm in the midline sagittal plane, but only at a few localized sites (primarily in the CSF near the edges of the thalamus and brainstem). In these regions, errors may be overestimated due to signal fluctuations caused by pulsatile CSF flow but have no effect on the ROI or region-label masks. Fig. 3 shows a plot of the cumulative number of selected points below a given level of error, demonstrating that approximately 90% of test points were perfectly aligned and that Fig. 2. Panel A shows a midline slice of the initial reference template, which is based on the average of one data set across all time points. An example of a single normalized midline sagittal slice is shown in Panel B, and Panel C shows the midline slice averaged over 24 experiments in eight volunteers. Fig. 3. Cumulative fraction of selected anatomical points as a function of the optimal calculated shift (i.e., error) to align midline sagittal and coronal planes with a reference volume. Errors are shown individually along each of the three axes, as well as the magnitude of error in the sagittal plane. The magnitudes of errors in the L/R direction in the coronal plane were too small to be detected with the assessment method used. essentially all (including those in the CSF) were accurate within a few millimeters. Fig. 4 shows a normalized region-label mask in midline sagittal and coronal planes. The same region-label mask is also shown after applying the inverse of the normalization transformation, allowing it to be mapped back onto the brainstem and spinal cord in the original image data.

Discussion
The effectiveness of the normalization procedure that we have developed for the spinal cord and brainstem is demonstrated by the correspondence of anatomical details across subjects. This can be directly observed by the fact that distinct anatomical features remain well defined, even when image data are averaged across a group of subjects with large anatomical variability (Fig. 2). Estimation of the normalization accuracy based on selected anatomical points shows that 90% of the points deviate from the corresponding point in the reference volume by less than 1 mm in the sagittal plane and that 93% are within 2 mm (Fig. 3). Normalization errors in the midcoronal plane were too small to be detected with this method and are essentially zero. Errors were larger in the A/P direction than in the L/R or R/C direction, corresponding to the direction of the largest variation across subjects (most likely owing to neck tilt and spinal cord curvature). The reference points with larger errors (N2 mm) in the sagittal plane tended to occur at localized regions around the thalamus and in the CSF-filled subarachnoid space and, thus, may reflect that the signal intensity was variable in these regions as a result of CSF motion.
While the image data are observed to be aligned to a good approximation, future studies may attempt to refine the proposed normalization algorithm by discriminating between gray and white matter regions. The alignment between spinal cord segments and vertebrae may also be accounted for in future iterations of the normalization, particularly when extended into the lumbar spinal cord where the alignment is more variable across people. Nonetheless, the normalization accuracy in the cervical spinal cord and brainstem is within the range that is acceptable for brain fMRI [19] and is adequate for subsequent application of more precise normalization methods such as ANIMAL (Automated Nonlinear Image Matching and Anatomical Labeling) [20].
The ability to apply an automatic ROI mask and define a standardized set of anatomical region-label maps is demonstrated for both normalized and original image data (Fig. 4). Alternative region-label maps can also be constructed to identify specific nuclei or regions, as needed. This is essential for automatic and/or objective assessments of where activity occurs in fMRI of the brainstem and spinal cord because fMRI data do not provide contrast to delineate brainstem nuclei or spinal cord gray matter. Reporting fMRI results, whether for research or for future clinical applications, requires referencing the areas of activity to a coordinate system or to anatomical features. Relying on visual assessments or poorly defined anatomy is subjective and may be prone to errors or bias. The results of the normalization procedure we have developed demonstrate that we can now objectively report the locations of areas of activity with an accuracy of a few millimeters. Therefore, it is now possible to create normative activity maps in the brainstem and spinal cord with a consistent format, allowing us to compare results within or between patient populations (with spinal cord injury or disease) and healthy volunteers.

Conclusions
The method that we have developed for 3D normalization and temporal alignment of spinal cord and brainstem fMRI data has been demonstrated to be effective. The accuracy of the normalization and the temporal alignment has been demonstrated to be within 2 mm or less in 93% of test points, with a mean accuracy of 0.3 mm. This advance permits unbiased comparisons of results across subjects, regardless of cross-subject variability in spinal cord size or curvature. Also, a standardized ROI mask has been created and voxelby-voxel region labels have been defined, either of which can be displayed in a normalized orientation or mapped back onto the original fMRI image data. Although further refinements may still be required to precisely align gray matter and white matter structures, and spinal cord segments at all levels of the cord, this development is a significant step toward making spinal fMRI a practical and reliable tool for research or clinical assessments of spinal cord function.