1. Introduction
Multiparametric MRI (mpMRI) of the prostate has emerged as a reliable non-invasive imaging modality for identifying clinically significant cancer, enabling selective sampling of high-risk regions with MRI-targeted biopsies, aiding in therapeutic planning and active surveillance [
1,
2]. A clinical mpMRI examination consists of the acquisition of T
2-weighted (T2w) images, DWI, and a gadolinium-based DCE series, deriving both anatomical and functional information from the prostate. It has been suggested that combining mpMRI with blood-sampled PSA-derived biomarkers could relieve uncertainty in the initial diagnosis [
3]. Radiomics in prostate cancer, focused on quantitative data from, e.g., MRI, could improve assessments of prostate cancer risk categories at initial diagnosis, allowing a personalized treatment approach [
4,
5].
Another quantitative MRI modality, proton magnetic resonance spectroscopic imaging (
1H-MRSI), has also been shown as useful because healthy and cancerous prostate tissue regions present different choline and citrate signal intensities [
6,
7,
8]. Ratios including these metabolites in combination with anatomical T2w MR imaging can identify and localize prostate cancer [
9]. Previous work in a multi-center setting has also shown a correlation between the metabolite ratios and prostate cancer aggressiveness classes [
10], indicating a role for
1H-MRSI in characterizing prostate cancer.
Clinical use of prostate spin–echo
1H-MRSI methods, however, is hindered by its limited robustness. Challenges in acquisition and post-processing techniques require special expertise [
11] and the relatively long acquisition time of 5 to 10 min might outweigh the additional value to routine multiparametric MRI [
12]. Next to these long acquisition times, the large voxel sizes of approximately 0.9 cm
3 to gain a sufficient metabolite signal in MRSI is regarded as a drawback. This low spatial resolution causes voxel bleeding to extend over larger parts of the prostate. Voxel bleeding is apparent from the side lobes of the spatial response function (SRF), describing the origin of the signal for a specific voxel of the MRSI matrix, which extends all over the field of view (FOV) of the MRSI grid [
13]. This spread of signal hinders precise localization of the signals of interest and causes the interference of other metabolites in the region of interest, such as peri-prostatic residual lipid signals extending into prostate tissue.
The side lobes of the SRF can be suppressed by applying a Hamming filter to the k-space data before spatial Fourier transformation (FT), decreasing far-reaching voxel bleeding [
14]. In combination with the weighted sampling of k-space, this is currently often used as an
1H-MRSI acquisition strategy [
15]. Applying this filter comes at the cost of an increase in the effective voxel size, seen as a widening of the SRF main lobe, leading to an even lower spatial resolution of the resulting metabolite maps [
14].
As an alternative to weighted sampling, an overdiscretized SRF target-based reconstruction method was presented by Kirchner et al. [
16]. Developed as an improved SENSE reconstruction method [
17] for single-slice two-dimensional (2D)
1H-MRSI brain data, the method was shown to be effective in suppressing near and far-reaching voxel bleed and correcting for intra-voxel coil sensitivity variations. In later work, a noise decorrelation was added to the method to correct for intra-voxel frequency shifts, improving spectral line shapes and increasing the signal-to-noise ratio (SNR) in the presence of the main magnetic (B0) field inhomogeneities on top of the voxel bleeding effects [
18]. These methods were applied to a single brain slice and limited to a spatial 2D field of view. The potential to extend the application to other regions within the human body studied with 3D-MRSI has not yet been explored.
In the prostate, a small organ in the center of the lower abdomen, the B0 field is mostly quite homogeneous, with occasional B0 field distortions at the posterior edge of the prostate resulting from the presence of rectal gas [
19]. Without large intra-voxel B0 inhomogeneities, the method of noise decorrelation with intra-voxel frequency shifts using B0 field maps is not expected to improve the spectra too much. We propose the use of random frequency shifts for the noise decorrelation of the spectra before the last step of applying a target function to reconstruct the final voxel sizes.
In this work we present the development of a three-dimensional overdiscretized and target-based MRSI reconstruction method (3D ODR method) to improve the spatial resolution of the 3D volume-selective MRSI of the prostate. After the initial overdiscretization of the MRSI grid, we used random frequency shifts for noise decorrelation before the last step of applying a target function. We compared this method with the conventional 3D MRSI with weighted k-space sampling and Hamming filtering of k-space at 3 Tesla. Without increasing the acquisition time from current 3D MRSI acquisition methods, the aim is to improve the localization of metabolite signals in the prostate without compromising on SNR.
2. Materials and Methods
2.1. Theory
The overdiscretized reconstruction method consists of three steps: overdiscretization of the MRSI grid, noise decorrelation between the subvoxels, and application of a target function to return to the final spatial resolution [
16,
18,
20]. For the current 3D volumetric application, multiple single-slice operations in the existing 2D approach were replaced by volumetric reconstructions, while parallel imaging was omitted, coil elements were combined by the scanning system, a spin–echo sequence is used, and the method of noise decorrelation is revised.
The reconstruction method starts with overdiscretizing the MRSI grid in k-space by zero-padding the original grid (of N voxels). This creates an intermediate grid with, after 3D FT, an increased interpolated spatial resolution. The size of the intermediate grid is determined by the choice of the overdiscretization factor (ζ). The same factor was applied to all three dimensions, resulting in a ζ
3-fold overdiscretized spatial grid. In target-driven overdiscretized reconstruction [
20], a minimization is performed of the cost function (Δ) for each voxel:
where F,
, and T are the reconstruction operator, the overdiscretized encoding matrix, and the target matrix, respectively. The analytical solution of the optimization problem is given by:
with
in which
reconstructs to the intermediate spatial grid of ζ
3N subvoxels. After applying
on the zero-padded k-space data s
OD,
the noise of the spectra’s neighboring subvoxels in the intermediate grid ρ
OD are strongly correlated. This correlation can be decreased, as the second step in the reconstruction, by a per-voxel frequency shift of the spectra (Δf(
r)) of the subvoxels at locations
r, resulting in a corrected intermediate grid:
The frequency shifts per voxel for the noise decorrelation were constructed with a map with random spatially distributed shift factors applied on the intermediate grid of the subvoxels. With the random spatial distribution of spectral shift factors, it is ensured that the noise between subvoxels is decorrelated, decreasing noise after weighted summation with a target function in the next step. With the number of spectral points to shift chosen within a small range (ranging from minus 2 to plus 2 spectral points), the widening of metabolite peaks is constrained.
As the third and final step, a 3D target function (T) is applied to the spatial intermediate grid after noise decorrelation to act similar to a spatial averaging operator:
This target function sums up spectra from the individual subvoxels in a weighted manner to return to the original grid resolution. This weighted spatial average of subvoxels, achieved by an SRF target with a 3D Gaussian shape that suppresses long-range signal bleeds, leads to the final effective voxel size of the reconstructed MRSI data, which is close to the original nominal voxel size.
A schematic representation of the 3D ODR method is shown in
Figure 1. The reconstruction method acts on two aspects: the overdiscretization and target function improve the SRF shape and corresponding signal localization while the noise decorrelation results in improved SNR. From previous work, the combination of an overdiscretization factor of 3 and target function width (σ) of 1.5 subvoxels was adopted, as this was shown to be optimal for SRF shape improvement [
16]. The SRF shape after ODR with these parameters has very low residual side lobes (
Figure 1).
With the chosen overdiscretization factor 3 and subsequent averaging with the target function of width 1.5, the resulting SRF has a width at a half height of 1.34 vs. 1.82 for the Hamming-filtered data (width relative to the nominal voxel width of the measurement) [
16]. These numbers can be extended from a 2D voxel diameter representation to a 3D voxel volume representation, using half this width in the formula for calculating the volume of a sphere (V = 4/3 × πr
3). An ODR dataset would represent a voxel volume of 1.26 times a nominal cubic voxel volume, as opposed to a voxel volume of 3.16 times the nominal volume for a conventional weighted dataset with Hamming filtering.
To illustrate the benefits of the improved SRF shape by our 3D ODR method we performed phantom measurements to analyze the effect of the reconstruction method on signal localization and measurements on a healthy volunteer to assess the feasibility in vivo for the application to 1H-MRSI of the prostate.
2.2. Acquisition
All measurements were performed on a Siemens 3T MR system (MAGNETOM Prisma-fit, Siemens Healthineers AG, Erlangen, Germany) using an 18-channel body phased-array coil and a 32-channel spine phased-array coil for signal reception (Siemens Healthineers AG, Erlangen, Germany). T2w images, used as reference and anatomical background, were acquired in three directions by a TSE sequence (TE = 101 ms, 2 averages; Tra: TR = 5660 ms, 0.3 × 0.3 × 3.0 mm; Sag: TR = 5590 ms, 0.6 × 0.6 × 3.0 mm; Cor: TR = 5000 ms, 0.6 × 0.6 × 3.0 mm). We acquired 3D
1H-MRSI with a GOIA semi-LASER
1H-MRSI sequence [
21] with scan parameters: TR = 750 ms, TE = 88 ms, spectral bandwidth = 2.4 kHz, and 1024 time samples. As the signal from citrate, the most abundant metabolite in healthy prostate, originates from a strongly coupled spin system, its spectral shape depends on sequence timing. For the semi-LASER sequence, the echo time of 88 ms and corresponding pulse timing results in an optimal signal shape of citrate in the prostate [
21].
The phantom used in this study is a cylindrical container filled with a sodium acetate solution with a cubic inner compartment with a creatine solution (
Figure 2). This box-shaped inner compartment has distinct straight edges and can therefore be used for localization analyses. The singlet of creatine at 3.04 ppm in the inner cube is the signal of interest and the
1H-MRSI is acquired with a matrix size of 11 × 11 × 11 voxels and FOV of 77 × 77 × 77 mm
3. We used two acquisition approaches: one intended for ODR, in which k-space is fully sampled in an acquisition time of 16:40 min with 1 average; the other, in which k-space is sampled in a Hamming-weighted manner [
14] with 8 k-space center averages to approximate a similar acquisition time (15:21 min). With these settings, both acquisitions yielded a similar acquisition time and matrix.
The acquisition strategy above does not result in equal voxel sizes for both datasets, as the width of the main lobe in the SRF is different for the weighted sampling and Hamming filtering versus the full k-space sampling with ODR. A second set of phantom measurements was aimed at acquiring similar final voxel sizes for both acquisition strategies within the same acquisition time. The mentioned SRF sizes, 1.34 vs. 1.82 times the radius of the spherical voxel for full k-space sampling with ODR and for weighted k-space sampling and Hamming filtering, respectively, were used to set up matrices resulting in equal voxel sizes. The first 1H-MRSI acquisition with a full sampling of k-space was performed with a matrix of 11 × 11 × 11 voxels, FOV of 88 × 88 × 88 mm3, and 2 averages in 28:52 min (voxel radius = 0.5 × 8 × 1.34 = 5.36 mm). The subsequent acquisition with weighted sampling of k-space was performed with 6 k-space center averages on a matrix of 15 × 15 × 15 voxels, with an FOV of 90 × 90 × 90 mm3 and acquisition time of 30:23 min (voxel radius = 0.5 × 6 × 1.82 = 5.46 mm).
Next to the two sets of phantom measurements, datasets from a healthy 45-year-old male volunteer, obtained in accordance with local ethics regulations, were included. No preparation to suppress peristaltic motion was performed. First acquisition with full k-space sampling (1 average) was acquired in 11:10 min with a matrix size of 11 × 9 × 9 voxels, FOV 88 × 72 × 72 mm3, MEGA water, and fat suppression. A second acquisition with weighted k-space sampling was acquired with 9 averages and an acquired FOV size of 77 × 63 × 63 mm3 (matrix size 11 × 9 × 9) in an acquisition time of 11:04 min. With these measurements, a similar acquisition time and matrix were achieved for both acquisitions.
Figure 2.
The phantom used in this study contains a creatine solution in the inner cube (spectrum from blue voxel, panel (B)) and acetate in the surrounding cylinder (spectrum from orange voxel, panel (C)). The matrix grid is indicated in yellow; the VOI is outlined in white (A). The colored lines in panel (A) show each plane in the other dimensions. The main peak of creatine (at 3.03 ppm) was set to represent the signal inside of the central cube and the integral of its fit was used to make creatine maps.
Figure 2.
The phantom used in this study contains a creatine solution in the inner cube (spectrum from blue voxel, panel (B)) and acetate in the surrounding cylinder (spectrum from orange voxel, panel (C)). The matrix grid is indicated in yellow; the VOI is outlined in white (A). The colored lines in panel (A) show each plane in the other dimensions. The main peak of creatine (at 3.03 ppm) was set to represent the signal inside of the central cube and the integral of its fit was used to make creatine maps.
2.3. Reconstruction
The weighted datasets were reconstructed by Hamming filtering of k-space, followed by zero-filling to a matrix of 16 × 16 × 16 voxels and spatial FT. This dataset functioned as a standard for comparison with the ODR pathway. The fully sampled datasets were zerofilled to 16 × 16 × 16 matrices and were processed by the 3D overdiscretized reconstruction method in Matlab (The Mathworks, Natick, MA). For the first step, an overdiscretization factor of 3 increased the MRSI grid of 16 × 16 × 16 voxels into an intermediate MRSI grid of 48 × 48 × 48 subvoxels. Noise decorrelation was subsequently performed per voxel on the high-resolution intermediate grid. Random spatial distributions of shift factors were simulated with the same resolution, e.g., 48 × 48 × 48, and frequency shifts were imposed on corresponding subvoxels. Shift factors were chosen in the range from −2 to 2 spectral points to ensure noise decorrelation while restraining spectral peak broadening. The noise-decorrelated intermediate spectroscopic image was transformed to the target resolution of 16 × 16 × 16 voxels by the SRF target function with σ = 1.5 subvoxels, as shown successfully in combination with the overdiscretization factor 3 in previous 2D work [
16].
2.4. Signal Analysis
The data processed with 3D ODR in Matlab and the originally acquired and postprocessed data (3D Hamming filter of zero-filled k-space) were imported in syngo.via (Siemens Healthineers AG, Erlangen, Germany) and all fitted with the same prior knowledge. For the phantom datasets, the signals of creatine (3.03 ± 0.05 ppm and 3.91 ± 0.05 ppm), lactate (1.34 ± 0.05 ppm), and acetate (1.9 ± 0.1 ppm) were fitted with a Gaussian function. The main peak of creatine (at 3.03 ppm) was chosen to represent the signal inside of the central inner cube and the integral of its fit was used to make metabolite maps. Creatine maps were normalized to the average integral value of the 8 voxels in the center of the cube. The creatine maps for the ODR data and Hamming-filtered data were compared both side to side as well as after subtracting one from the other in a difference map for a quantitative analysis of differences in voxel bleed. SNR was defined as the peak integral of creatine per voxel divided by the noise, stated as the standard deviation of the signal in an area of the spectrum not containing metabolites. To compare SNR values between the datasets with different voxel sizes, only voxels completely within the inner cube were chosen. These center 8 voxels overlapped with each other; the SNR values were therefore not independent. As the SNR is determined by both the SRF shape and the frequency shift, the separate contributions will be further analyzed by rendering an additional dataset for each experiment of full k-space sampling by applying ODR without the per-voxel frequency shift.
For the in vivo datasets, signals from citrate and residual lipids were fitted with appropriate line shapes at 2.6 ± 0.1 ppm and 2.3 ± 0.2 ppm, respectively. For the combination of total choline (3.21 ppm), polyamines (spermine at 3.12 ppm), and creatine (3.03 ppm) signals, the individual metabolite peaks were fitted as separate peaks. Metabolite maps of citrate and lipids were created and used for a qualitative analysis of signal localization. SNR values were extracted from the integral of the citrate peak over the noise value, where noise was again defined as the standard deviation of the signal in an area of the spectrum not containing metabolites. To compare spectral quality between ODR with and without frequency shift and Hamming-filtered data, SNR values from all datasets were compared in a voxel-wise manner, with a correction to account for the difference in relative voxel sizes. A selection of voxels from the same spatial location in all datasets, which were completely within the prostate with minimal overlap with each other, were chosen for SNR comparison. Localization of metabolite signals in both datasets was qualitatively analyzed by comparing citrate and lipids maps.
4. Discussion
In this work, we developed a 3D overdiscretized reconstruction method for semi-LASER 3D 1H-MRSI and applied it to phantom and in vivo prostate measurements. With the 3D ODR method, it is possible to suppress far-reaching voxel bleed, with less increase in the effective voxel size as with weighted k-space sampling and Hamming filtering. In addition to the superior signal localization, a moderate gain in SNR was achieved.
Spatial resolution is an important aspect in MRSI, as a relatively large voxel size is needed to acquire signals from metabolites instead of water. This gives rise to partial volume effects, with one voxel representing a signal from different tissues, as is visible in the SRF, which hinders detecting or characterizing small lesions. With Hamming filtering, spatial resolution is traded for the suppression of far-reaching voxel bleed, with a widened main lobe of the SRF and decreased side lobes. The 3D ODR method increases the true spatial resolution of metabolite maps in comparison with Hamming filtering, matching underlying anatomical structures as found on clinical prostate MR imaging, and may therefore improve the metabolic characterization of small prostate cancer lesions.
For the second step of the ODR method, noise decorrelation of the overdiscretized 3D MRSI grid, two strategies can be used. The initial overdiscretization into an intermediate grid is created in real space by the zero-filling of the acquired k-space. Overlap of the respective SRFs is the major source of noise correlation in neighboring subvoxels of the intermediate grid. Imposing a small random frequency shift on the subvoxels before the third step of applying a target function to return to the final spatial resolution leads to noise decorrelation. If noise decorrelation is based on measured B0 maps, the frequency alignment of subvoxels can improve linewidths. However, with noise decorrelation based on random spectral shifts, instead of using a B0 map, spectral line widths can increase [
18]. In phantom measurements with well-resolved and very narrow spectral lines, this spectral line broadening is visible. Choosing only few random spectral points to shift relative to the linewidth of the metabolites of interest prevents strong line broadening. Linewidths were larger in measurements of the prostate than in the phantom, not only because of tissue susceptibility differences but possibly also due to small body movements [
22]. Noise decorrelation in the intermediate grid was based on random spectral point shifts between −2 and +2 points. The resulting spectra from the volunteer measurements showed some peak broadening. However, its impact on metabolite analysis is limited with the usual processing of quantifying the total spectral integral from 3.21 to 3.04 ppm with the sum of three spectral shapes of choline, spermine, and creatine, and with only a loss of the splitting of the two center peaks of the strongly coupled spectral shape of citrate.
Application of the ODR method, both to the phantom and volunteer measurements, showed an improved localization of signals of interest and decreased lipid contamination within the prostate of the volunteer. Metabolite localization and SNR could be directly compared in the phantom between weighted k-space sampling with Hamming filtering and full k-space sampling with ODR, as the metabolite concentration and localization within the inner cube of the phantom was homogeneous and demarcated. The resulting difference maps in all three dimensions visualized the improved signal localization as a result of the improved 3D SRF shape when measuring with the same acquisition matrix in the same acquisition time. For the volunteer measurements, overlaying both datasets on anatomical T2w images provided the location of the peripheral zone or whole prostate gland as boundaries for estimated signal localization. However, comparing the signal intensities of voxels with different sizes remains challenging, as tissue metabolite signal intensity is not necessarily homogeneous, with varying metabolite concentrations across the whole prostate. Some studies have shown the use of anatomical boundaries to constrain metabolite signals within these boundaries, with the risk of biases in quantifying differences in metabolite signals within this anatomical region [
23]. Apart from voxel definition by the two ways of MRSI data processing, accurate localization of metabolite signal within the whole volume of interest is also due to the use of a semi-LASER pulse sequence in the
1H-MRSI acquisition of the phantom and volunteer. This volume selection already provides more accurate localization over a PRESS sequence for the whole VOI [
21].
The effect of ODR on the observed SNR in the two separate phantom measurements needs clarification. An increase in SNR of 1.87 was achieved for the fully sampled and ODR dataset compared with the weighted k-space sampling and Hamming filtering data in the first measurement. This measurement involved different voxel sizes and a factor to correct for this choice. In the second phantom measurement we aimed at achieving a similar voxel volume for both the ODR dataset and the weighted Hamming-filtered dataset. Matching final voxel sizes was successful, as could be seen from the localization difference image between the two acquisitions (
Figure 3F). The ODR method resulted in a SNR increase of 45%. In earlier work, Nassirpour et al. [
20] found an average 23% SNR gain in ODR data with SENSE compared with reconstructing with SENSE alone (with ζ = 4 at a field strength of 9.4T). Kirchner et al. [
18] found an increase in SNR between the ODR data and the Hamming-filtered data with a factor of 3.7 (with ζ = 5 at a field strength of 7T). There are several differences between those previous studies and the current one, and next to field strength and B0 homogeneity in the FOV, the spatial resolution of the MRSI datasets, the exact definition of SNR (peak or integral), and the upgrading to 3D can influence the SNR gain.
For the SNR comparison in the first phantom measurements, we calculated 3D voxel sizes extrapolated from simulated 2D SRF values, assuming spherical voxel sizes with a diameter of FWHM of the SRF. Moving from a 1D sinc-like, to a 2D circular shape, to a 3D spherical SRF facilitates calculations for effective voxel volume in 3D from the simulated numbers for 1D and 2D. The effective final voxel volumes were subsequently used as the correction factor for SNR comparison. SRF line shape differences can cause differences in assumed voxel volumes between the two phantom measurements, resulting in different values for the correction factor for both signal and noise in the SNR of different voxel volumes.
In the datasets from the volunteer examination, it is more difficult to assess SNR gain. As the two datasets are of different spatial resolutions, one cannot choose the same position and voxel size in the two acquisitions; therefore, local tissue concentration differences can obscure a fair SNR comparison (even with correction for voxel size). SNR averaged over voxels completely within the peripheral zone was the best option here to compare both datasets, with full k-space sampling with ODR exceeding the SNR of the metabolites of interest over weighted k-space sampling with Hamming filtering.
Subsequent testing of the ODR method in patients with prostate cancer can establish the effect of improved signal localization and SNR for tumor characterization. With applying ODR on a fully sampled k-space dataset, the spatial resolution acquired in our in vivo examination required an acquisition time of 11 min. With expected differences in prostate sizes and larger prostate VOIs in older subjects, the acquisition time for MRSI might become too long. A solution would be to combine complete k-space sampling and ODR with accelerated read-outs, e.g., EPI or spiral read-out methods [
24].
With an accelerated acquisition and automated post-processing workflow, resulting high-quality spectra could be of additional value in the clinical management of prostate cancer. Especially in repeated examinations per patient, as in active surveillance, we believe acquiring quantitative MRSI data enables direct comparison between scans and maximizes non-invasively acquired information of the prostate.