Hitchhiker's Guide to Voxel Segmentation for Partial Volume Correction of In Vivo Magnetic Resonance Spectroscopy.

Partial volume effects have the potential to cause inaccuracies when quantifying metabolites using proton magnetic resonance spectroscopy (MRS). In order to correct for cerebrospinal fluid content, a spectroscopic voxel needs to be segmented according to different tissue contents. This article aims to detail how automated partial volume segmentation can be undertaken and provides a software framework for researchers to develop their own tools. While many studies have detailed the impact of partial volume correction on proton magnetic resonance spectroscopy quantification, there is a paucity of literature explaining how voxel segmentation can be achieved using freely available neuroimaging packages.

Introduction

Proton magnetic resonance spectroscopy (MRS) recorded from the brain allows for noninvasive quantification of tissue metabolites, using a conventional MRI scanner. In the brain, MRS has been applied to different diseases such as malignancy, multiple sclerosis (MS), cognitive impairment, and stroke.1 MRS provides information on tissue metabolism, in both health and disease, and has the potential to provide biomarkers for pathologies not traditionally appreciated on conventional imaging.2–4

To obtain single-voxel MRS datum, a three-dimensional volume known as a voxel is placed in a region of interest in the brain. A voxel is usually made up of the following three primary quantities: white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF). The spectroscopy voxel should not be confused with the imaging voxel, which represents a much smaller volume, typically 1 mm3 made out of WM, GM, or CSF, compared with at least 1 cm3 for single-voxel spectroscopy (SVS). The area under the metabolite spectral resonances that are obtained from within spectroscopy voxel is directly proportional to the in vivo concentration of metabolites.5 However, the metabolite signal intensity is affected by many factors such as radio-frequency (RF) coil loading and excitation homogeneity. To correct for these RF effects, metabolite peak intensities need to be normalized against an internal or external standard; this is referred to as relative quantification. Typically, metabolite concentration is expressed as a ratio to an internal reference such as total creatine (tCr), total N-acetylaspartate (tNAA), or water. However, if a change is noted in the concentration ratio, it can be difficult to ascertain if the change is due to the denominator, numerator, or both.

An alternative to relative quantification is absolute metabolite quantification. The most common method (water signal referencing) for absolute quantification is reliant upon tissue water concentration as an internal reference.6 The concentration of water varies according to the tissue compartments contained within the spectroscopic voxel. Variations in tissue and CSF contents within the voxel affect quantification in two main ways: the relative water content varies according to the tissue type (ie, CSF, which is mostly made out of water and has a greater concentration of pure water compared with GM and WM)7 and the majority of metabolites are contained within the GM and WM (partial volume effects). Therefore, partial volume effects should be considered for accurate quantification. If the water concentration is not corrected prior to metabolite quantification, metabolite concentrations will be underestimated.8 Despite partial volume correction being a recommended step for absolute quantification of metabolites,5 there is a paucity of literature describing precisely how voxel segmentation can be achieved from a practical standpoint. The following work was developed on a Prisma scanner (3 T; Siemens), and as a result, modifications might be required for general use. As part of this work, a Matlab script that performs all the steps detailed below in an automated fashion was developed (the script is available with the corresponding author).

Multiparametric MRI can be useful in increasing the sensitivity and specificity of a biomarker or to confirm MRS findings. Therefore, it is useful to have an understanding of how nonspectroscopy parameters (diffusion weighted imaging, functional MRI, and structural MRI) can be determined within the MRS voxel. The basic steps required to determine alternate MR metrics from within an MRS voxel will be outlined.

Co-registration of Voxel Mask

A binary mask must first be created, which represents the 3D SVS voxel. A binary mask is an image where each voxel contains either 1 or 0; in this case, the SVS voxel would be represented by a value of 1 and all other voxels 0. There are multiple ways to create a binary mask; one example, which was used successfully in this work, is the “mask()” function within Matlab9 and the FreeSurfer (Fischl)10 application “mri_volsynth.”

Prior to performing voxel segmentation, a binary mask representing the spectroscopy voxel must be co-registered to the imaging coordinates (i, j, k). The pixels within the voxel mask are transformed using an affine transformation (A), to generate scanner coordinates (x, y, z) according to Equation (1):

We can expand Equation (1) as shown below in Equation (2):

The rotation (r) and translation matrices (t) can be combined to form a transformation matrix (M):

Equation (3) can be used to transform structural image from scanner coordinates (x, y, z) to voxel coordinates (i, j, k) with the following transformation matrix:

Additionally, we can transform the scanner coordinates to MRS voxel coordinates (i′, j′, k′) using the following transformation matrix:

Finally, Equations (4) and (5) can be combined to create a translation matrix from anatomical space to SVS space according to Equation (6):

The transformation matrix is calculated by taking the product of the inverse anatomical image (MMRI)−1 transformation matrix and the SVS voxel (MMRS) transformation matrix. The (M)−1 describes the realignment of the voxel coordinates to scanner coordinates and (MMRS) describes the transform of the scanner coordinates to the MRS voxel.

The transformation matrix for the imaging datum (MMRI) can be obtained using the SPM command “V = spm_vol (nii_file),” which will create a structural array containing image volume datum.11 The “.mat” element of the “V = spm_vol (nii_file)” array contains the 4 × 4 affine transformation that maps voxel coordinates to scanner coordinates for the anatomical image. The transformation matrix for the SVS voxel (MMRS) is more difficult to construct, and the input datum is taken from the scanner exported “.rda” file. Table 1 outlines the variables required for the “.rda” file and their description.

Table 1

A description of the variables contained within the Siemens .rda spectroscopy file and how they relate to the voxel dimensions, translation, and rotational matrices.

VARIABLE NAME WITHIN.rda FILE	DESCRIPTION
VOIPhaseFOV, VOIReadoutFOV and VOIThickness	Voxel dimensions
VOIPositionSag, VOIPositionCor, VOIPositionTra	SVS voxel translation
ColumnVector [0] RowVector [0] ColumnVector [0] × RowVector [0]	Rotational matrix
ColumnVector [1] RowVector [1] ColumnVector [1] × RowVector [1]
ColumnVector [2] RowVector [2] ColumnVector [2] × RowVector [2]

The above steps were performed using a Matlab function created by Harris et al;12 however, it is possible to write a custom Matlab programme using the steps mentioned above. It is highly recommended to save screen shots of the voxel location on the MRI scanner for each SVS location, which can be achieved by using “alt-print screen” key combination and then pasting the imaging into Microsoft Paint on the MRI scanner computer. These reference images can then be used to ensure that the voxel mask created using the steps above has been correctly transposed into scanner coordinates.

Figure 1 shows an example of an SVS mask overlaid onto a T1-magnetization-prepared rapid gradient echo (MPRAGE) image after the SVS voxel has been transposed to the scanner space; for comparison, the SVS voxel prescription image at the time of MRI scan is also shown.

Figure 1

(A) SVS voxel prescription image taken at the time of scanning. (B) Reconstructed SVS voxel displayed as a mask overlying the T1-MPRAGE image. T1-MPRAGE image was acquired at 3 T using a 64-channel head coil: (TR/TE/TI = 2000/3.5/1100 ms, flip angle = 7°, field of view = 256 × 256 mm, image voxel size 1 × 1 × 1 mm3).

Segmentation of Voxel Prescription Image

To determine the quantities of tissue subtypes within an SVS voxel, a good quality structural image must be obtained at the time of scanning. A 3D T1-weighted MPRAGE sequence13 allows for identification of anatomical structures at the time of voxel prescription and voxel segmentation at the time of postprocessing. The following are typical parameters at 3 T for a 3D T1-weighted MPRAGE sequence: (TR/TE/TI = 2000/3.5/1100 ms, flip angle = 7°, field of view = 256 × 256 mm, image voxel size 1 × 1 × 1 mm3, IPAT 2, orientation = sagittal, acquisition time 4:48 minutes). Even though some other brain segmentation packages might require T1- and T2-weighted images, the above MPRAGE will suffice for the following procedure.

There are four major open-source neuroimaging packages that allow for segmentation of a 3D brain image into GM, WM, and CSF: FSL,14 FreeSurfer (Fischl),10 SPM,11 and AFNI.15 Unfortunately, FreeSurfer does not segment sulcal CSF,16 which could result in inaccuracies in partial volume correction. FSL FAST, which will be referred to as FAST,17 was chosen for this work; however, SPM or AFNI could also be used.

FAST segments tissue according to a hidden Markov random field model and an associated expectation–maximization algorithm, which has been shown to be robust and reliable when compared with finite mixture-based models that can be sensitive to noise.17 Prior to segmentation with FAST, the structural image should first be brain extracted using the brain extraction tool (BET)18 within FAST. After segmentation, FAST will output partial volume maps for GM, WM, and CSF as shown in Figure 2.

Figure 2

The results of partial volume segmentation using FAST (parameters described in text), (A) CSF mask, (B) WM mask, and (C) GM mask.

Determination of Tissue Fractions within SVS Voxel

To determine the fractions of WM, GM, and CSF within each SVS voxel, the voxel mask that has been registered to the anatomical space needs to be overlaid with the partial volume maps created using FAST (Fig. 2). Both FSL and FreeSurfer have utilities that will achieve this goal, in FAST the utility is “ fslstats.” Within FAST stats, the mean of the partial volume map for GM, WM, and CSF is determined, using the SVS voxel mask. The commands required to achieve this step are shown below in “Example commands 1—Quantifying tissue classes within an SVS voxel.”

Absolute Quantification with Partial Volume Correction—1D Spectroscopy

To quantify concentrations of metabolites, peak area must first be referenced to a known value. Referencing corrects arbitrary amplification of signal for errors associated with variable RF coil gain and allows comparisons to be made between peak intensities of different metabolites. There are several methods available for metabolite quantification such as internal metabolite referencing (typically to tCr or tNAA), internal water signal referencing, and external reference solutions.5 The most popular method for absolute quantification of metabolites using MRS is internal water signal referencing, which involves acquiring an unsuppressed water spectrum after the 1D MRS acquisition. It is essential that the unsuppressed water spectrum should be acquired in a manner that is identical to the water-suppressed spectrum, ie, same voxel location.5

When using water signal referencing, it is assumed that the fully relaxed water signal will be proportional to the number of water moles in the voxel.19 Additionally, it is assumed that CSF has no significant contribution to metabolite concentrations. Thus, the absolute metabolite concentration, M, is given by: where IMet is the observed metabolite signal, I(GM&WM) is water signal arising from GM and WM, #H is the number of protons contributing to the metabolite peak being quantified (eg, it is 3 for NAA peak at 2.02 ppm), 2 is the number of protons contributing to the water signal and [H2O] is the molar concentration of pure water (55.6 mol/L).

The I term can be derived the from the observed water intensity (I(obs)), accounting for the fact that the water signal arises from multiple compartments and is not fully relaxed: where f are the water fractions (note not the volume fractions) contained within each tissue class (GM, WM, and CSF), refers to the fully relaxed water signal, the relaxation attenuation factors ( ) are found using Equation (9) and accounts for the fact water relaxes differently depending on which tissue class is it is contained within: where TE and TR are the echo and repetition times, respectively, T1 and T2 are the T1 and T2 relaxation times of water in tissue compartment X.

Values for T1 and T2 relaxation times of water at 3 T can be found in the literature20,21 and are shown in Table 2.

Table 2

In vivo T1 and T2 water relaxation times (ms) at 3 T for relaxation correction using Equation (9).

	GM	WM	CSF
T1	1470	1600	3000
T2	110	74	200

The I(GM&WM) is found by rearranging Equation (8) for and subtracting the CSF fraction (fCSF):

Finally, Equations (7) and (10) can be combined to yield:

The fractions of water in different tissue compartments are represented by fGM, fWM and fCSF, which can be expressed in terms of tissue segmentation fractions: where, c is the water content in GM (0.82), WM (0.73), and CSF (0.98)22 and ν is the tissue volume fraction that is output as a result of SVS voxel segmentation. Ernst et al first showed that tissue compartments could be estimated using bi-exponential fitting of the T2 decay of brain water signal.7 However, for many studies, estimating the brain tissue water and CSF using relaxation experiments is not practical due to time limitations. As neuroimaging software packages have improved, the standard method of estimating the tissue compartments has become structural image segmentation.19,23,24

A common tool for metabolite quantification of 1D spectra is LCModel, where prior knowledge is used to fit a basis set of experimental or simulated spectra to the acquired in vivo spectrum.25 This technique is useful as it allows for the estimation of overlapping resonances such as glutamate and glutamine.5 Partial volume correction can be applied using LCModel by adjusting the “WCONC” term according to the equation8,26:

If the above approach is used, “WCONC” variable needs to be calculated for each SVS voxel prior to LCModel analysis and input at the time of fitting.

Alternatively, total water concentration in the voxel can be ignored and a correction for CSF can be applied according to27: where C is the corrected metabolite concentration, C0 is the LCModel output and νCSF is the volume fraction of CSF contained within the SVS voxel. This method only considers the partial volume effect of CSF.

Partial Volume Correction—2D Spectroscopy

Two-dimensional (2D) MRS allows composite or overlapping resonances from one-dimensional (1D) spectra to be separated out. In conventional 1D spectroscopy, intensity (y-axis) is plotted against frequency (x-axis), whereas in 2D spectroscopy, intensity is plotted against two frequency variables.28 2D spectroscopy was first introduced by Jeener29 as a two-RF pulse sequence known as correlated spectroscopy (COSY).30,31 The introduction of 2D in vivo spectroscopy has allowed researchers to make unambiguous assignments that could not have been made using 1D spectroscopy and has also increased the clinical utility MRS.32–34

Unfortunately, 2D MRS quantification methods are not as mature as those available for 1D MRS. Typically, 2D datum is quantified by integration after manual peak picking as described here;35 however, prior knowledge-based fitting algorithms have started to emerge.36 The most common peak quantification method used when analyzing 2D spectroscopy is the internal endogenous marker method. When using this technique, a peak ratio is calculated using an internal reference (for example, diagonal singlet of tCr or tNAA), the metabolite concentration is then determined by multiplying the peak ratio by the in vivo concentration of the internal reference.5 However, most commonly, the final step is not carried out and just the peak ratio is determined.

It is possible to reference peak volumes to the internal water signal; the main advantage is that the water signal is less affected by physiological perturbations than other metabolites. However, this technique still needs to be used cautiously as several studies have found that brain water can change in conditions such as MS,37 neoplasia,38 and hydrocephalus.5,38

Thus, taking partial volume correction into account, the metabolite ratios are determined according to: where fGM, fWM, and fCSF can be determined using Equation (12).

Quantification of Other MRI Metrics

Being able to quantify other MRI metrics within an SVS voxel can be useful for carrying out correlation or regression analysis. For example, as a part of a MS study, it may be of interest to know the total volume of T2 hyperintense lesions within the SVS voxel. Alternatively, researchers may like to correlate diffusion metrics such as fractional anisotropy (FA, a measure of the degree of diffusion anisotropy) or mean diffusivity (MD, a measure of apparent diffusion).39

We will assume that diffusion datum has been corrected for eddy current distortions and a tensor model has been fitted, allowing FA and MD to be calculated. Prior to extracting FA and MD values within the SVS voxel, the diffusion datum must be co-registered to the structural image that was used to prescribe the SVS voxel location. This can be achieved using the FAST tool FLIRT,40,41 specifically the “epi_reg” script. This script registers diffusion EPI images to structural images, such as T1-weighted images. In the example commands shown below, the b = 0 seconds/mm2 images are registered against the structural images, as they provide better tissue contrast, which is important for good quality image registration. Additionally, in this example command, “epi_reg” will perform EPI distortion correction at the same time as EPI registration. After registering the diffusion tensor images (DTI) to the structural space, the mean values for MD or FA can be determined using the “fslstats” command. A summary of this process is shown below in the “Example commands 2” box.

This same process of co-registration, transformation, and extraction can be used for any other imaging metric where datum is prescribed to an imaging voxel, eg, blood oxygenation level-dependent imaging. However, if the metrics of interest, eg, WM lesions in a MS study, are already in the structural space, then the registration step can be skipped.

Example Registration and Segmentation

A single-voxel 1D spectrum was acquired from the hippocampus using a voxel size 3 cm × 1.5 cm × 1.5 cm = 6.75 cm3, from a patient who had previously been diagnosed with relapsing and remitting multiple sclerosis. MS is a neurodegenerative condition characterized by autoimmune-mediated demyelination. At the time of scanning (Prisma; Siemens), the following structural images were acquired using a 64-channel head coil: 3D T1-weighted MPRAGE sequence (TR/TE/TI = 2000/3.5/1100 ms, flip angle = 7°, field of view = 256 × 256 mm, image voxel size 1 × 1 × 1 mm3, IPAT 2, acquisition time 4:48 minutes); T2 fluid-attenuated inversion recovery (FLAIR; TR/TE/TI = 5000/386/1800 ms, echo train duration: 858 ms, field of view = 256 × 256 mm2, with spatial resolution of 1 × 1 × 1 mm3, IPAT 3, acquisition time 4:12 minutes).

The SVS voxel mask was created according to the steps outlined above in the “co-registration of voxel mask section.” The results of this step are shown in Figure 1 and can be compared with the original prescription voxel image that was saved at the time of scanning.

On T2 FLAIR, MS plaques appear hyperintense, and thus, it is important to know if there is a lesion within the SVS voxel, as the spectrum from these areas will be different to normal appearing brain.1 Figure 3 shows an example of a hippocampal SVS voxel that contains WM lesions. T2 FLAIR lesions were segmented using the lesion growth algorithm as implemented in the LST toolbox version 2.0.6 for SPM.42 The algorithm first segments the T1 images into the three main tissue classes (CSF, GM, and WM). This information is then combined with the co-registered FLAIR intensities in order to calculate lesion belief maps. By thresholding these maps with a prechosen initial threshold (κ), an initial binary lesion map was obtained, which was subsequently grown along voxels that appear hyperintense in the FLAIR image. The initial κ threshold was selected by iterating κ and performing a visual inspection, for this data a κ value of 0.1 was selected. A binary lesion mask was created for each participant using a threshold of 0.5. Hypointense lesions on the T1-MPRAGE were determined using the binary lesion mask and filled with intensities similar to voxels not contained within a lesion. Lesion filling was performed using the LST toolbox42 to improve volume measurements and prevent errors in partial volume segmentation.43

Figure 3

T1-MPRAGE image of patient with MS, the hippocampal SVS voxel is overlaid in red, and the binary lesion map that was constructed using a FLAIR image is shown in yellow. Notice that there are WM lesions contained within the hippocampal SVS voxel. T1-MPRAGE image was acquired at 3 T using a 64-channel head coil: (TR/TE/TI = 2000/3.5/1100 ms, flip angle = 7°, field of view = 256 × 256 mm, image voxel size 1 × 1 × 1 mm3).

Tissue classes were segmented from the lesion-filled T1-MPRAGE using FAST17 as described above. The tissue class volume fractions within the hippocampal voxel were calculated according to the steps in the “Example commands 1” box above. Finally, the volume fractions and T2 FLAIR lesions were quantified within the hippocampal voxel; the results of this step are shown in Table 3.

Table 3

Fractions of CSF, GM, WM, and WM lesion volume contained within the example SVS voxel.

	CSF	GM	WM	WM LESION (mm3)
Hippocampus	0.029	0.421	0.550	117

Conclusion

Partial volume effects have the potential to introduce error when quantifying metabolites using proton magnetic resonance spectroscopy (MRS). Herein, we have provided a comprehensive outline of how to undertake partial volume correction using freely available neuroimaging packages and how these corrections can be applied to 1D and 2D spectroscopy (L-COSY). Additionally, a guide to extracting other MRI metrics of interest (such as FA or MD) from the SVS voxel volume was provided.

Extract the brain from the

$ t1image.nii extractedbrain.nii.gz -f 0.3 -g 0

Perform tissue segmentation

$ fast -t 1 -n 3 -H 0.1 -I 4 -l 20.0 -o /fast_output_folder/extractedbrain.nii.gz

Determine the grey, white and CSF content of the voxel

$ fslstats -t /fast_output/fast_output_pve_0 -k SVS_mask.nii –m

Note: ‘fslstats’ (step 3 above) needs to be executed for each tissue type, these masks are located in the fast_output folder (see step 3 above) and are named pve_0, pve_1 and pve_2 (shown in Fig. 2 is this order) assuming three tissue types have been used.

Co-register diffusion datum with SVS prescription structural image

$ epi_reg –-epi=my _hifi_b0 –t1=t1image.nii – t1brain=extractedbrain.nii.gz –out=dwi2struct –fmap my _fieldmap_rads –fmapmag=my _fieldmap_mag – fmapmagbrain my_fieldmap_mag–echospacing=0.000345 – pedir=-y –wmseg/fast_output/fast_output_pve_2

Transform the DTI datum into the structural space

$ applywarp –in=FA_map.nii.gz –ref=t1image.nii –out=FA_2_ struct –warp=dwi2struct_warp

Determine the mean value for FA or MD within the SVS voxel

$ fslstats -t FA_2_struct.nii.gz -k svs_mask.nii -m