Distortion‐free 3D diffusion imaging of the prostate using a multishot diffusion‐prepared phase‐cycled acquisition and dictionary matching

To achieve three‐dimensional (3D) distortion‐free apparent diffusion coefficient (ADC) maps for prostate imaging using a multishot diffusion prepared‐gradient echo (msDP‐GRE) sequence and ADC dictionary matching.


| INTRODUCTION
Multiparametric MRI (mpMRI) is now recommended for the detection and risk stratification of prostate cancer, as demonstrated by recent guidelines, 1,2 based on data from multicenter studies. [3][4][5][6] Diffusion-weighted imaging (DWI) is a key component of prostate mpMRI for its ability to detect cancerous lesions. 7,8 Current clinical DWI protocols are based on two-dimensional (2D) single-shot diffusion weighted-echo planar imaging (ssDW-EPI), which is fast to overcome the issue of motion corruption, but of limited spatial resolution and prone to geometric distortion caused by off-resonance susceptibility. This limits the diagnostic performance of diffusion MRI and prohibits its use in emerging MR-guided precision treatments, eg, in radiotherapy planning with its high requirement for 3D geometric accuracy. Thus, there is a clear need to improve the implementation of this key cancer imaging contrast.
Multishot sequences can achieve high and isotropic resolution and have been successfully used for DW-EPI aiming to reduce geometric distortion, [9][10][11][12] but the acquisition often takes too long and the images may still be prone to distortion. Multishot sequences can also be combined with diffusion preparation (DP) modules (as opposed to DW), where the diffusion-sensitized transverse magnetization is stored into longitudinal magnetization by a 90° tip-up pulse. DP sequences have the advantage that diffusion encoding and readout are separated, and can thus be independently optimized for resolution, distortion, and signal-to-noise ratio (SNR), and can be combined with any readout scheme: turbo spin echo (TSE), 13,14 balanced steady-state free-precession (bSSFP), 15,16 and gradient echo (GRE). 17,18 In multishot diffusion-prepared (msDP) sequences, there are two main challenges that cause inaccurate apparent diffusion coefficient (ADC) quantification.

| Challenge 1: magnitude errors
The first major issue is accurate diffusion quantitation in the presence of additional (unknown) signal phase, which results in signal magnitude loss as magnetization is tipped up at the end of the DP module, thus confounding the desired diffusion weighting (challenge 1). These magnitude inconsistencies can be caused by any unwanted motion (patient movement, cardiac pulsation, table vibration) occurring while the large diffusion-encoding gradients are applied. Typical sources of motion in prostate imaging are bowel motion, and bladder and rectal filling. These intershot motion-induced magnitude errors (MiMEs) are (1) spatially dependent and (2) shot dependent, leading to magnitude inconsistencies when multiple corrupted shots are combined together in k-space. Several methods have been introduced to account for intershot MiME: motion-compensated DP modules to null the firstorder or second-order gradient moments, 19,20 and phase cycling of the DP tip-up pulse to recover the full transverse magnetization in two separate acquisitions, sampling both the x and y components of the transverse magnetization 17 ; however, this approach doubles the scan time and might not work if motion between the subsequent acquisitions is considerably different. Additional dephasing gradients before the DP 90° tip-up pulse (magnitude stabilizers) have been proposed to create a uniform-phase distribution throughout the shots and stabilize the signal magnitude. 21 As the magnitude stabilizers convert the magnitude inconsistency back to phase inconsistency, these can be used in concert with navigated intershot phase-correction methods. 13,14,16 However, the magnitude-stabilizer approach comes with a 50% signalloss penalty, thus further limiting the already critical SNR efficiency in diffusion imaging.
Another source of magnitude errors (MEs) is the presence of eddy currents arising from the fast switching of the strong diffusion gradients. These eddy current-induced magnitude errors (EiMEs) are (1) spatially dependent, (2) constant throughout the multishot acquisition, and (3) proportional to the strength of the diffusion gradients. The eddy current problem has been previously accounted for by employing a pair of bipolar rather than monopolar gradients to compensate currents with opposite signs, [22][23][24][25] by performing phase cycling of the DP tip-up pulse, 17,25 or by adding magnitude stabilizers gradients. 20

| Challenge 2: T 1 effects
Another quantitation challenge in rapid msDP sequences is that the acquired signal is a mixture of DP and T 1 -recovered magnetization (challenge 2). This is caused by several effects: (1) the shot length (pulse repetition time [TR]) being shorter than the underlying tissue T 1 , which prohibits full recovery of the longitudinal magnetization between shots; (2) the longitudinal magnetization regrowth during readout; (3) the unprepared (not DP) magnetization caused by an imperfect flip angle (FA) of the DP 90° tip-down and tip-up pulse, such that some of the magnetization will not experience the diffusion | 1443 ROCCIA et Al. gradients and will start regrowing. These T 1 effects make a simple monoexponential fitting unsuitable for accurate ADC quantification. 17,24 Possible solutions that have been suggested to account for these T 1 -related effects are: to increase TR or add a delay time between the shots to allow for full magnetization recovery, at the cost of increased acquisition time 17,24,25 ; to use a center-out trajectory, as the first readout segments are less T 1 -contaminated than later segments 18,19,24 ; to use phase cycling 25 ; or to use magnitude stabilizers, 13 which do not refocus unprepared magnetization during data acquisition.
The aim of this study was to address the aforementioned challenges to obtain nondistorted 3D ADC maps of the prostate. We proposed a 3D msDP sequence with a high number of shots, using a GRE readout and a custom Cartesian, centric 3D k-space trajectory with center oversampling. To address challenge 1 (EiME and MiME), a k-space center averaging and phase cycling strategy were proposed. The centeroversampling property of the trajectory was exploited for motion-averaging across many shots, while undersampling the periphery of k-space combined with total variation sensitivity encoding (tvSENSE 26 ) reconstruction allowed shortening of the acquisition time. To address challenge 2 (T 1 effects), we used a centric trajectory and ADC dictionary matching based on extended-phase-graph (EPG 27 ) signal simulations.

| METHODS
All data were acquired on a 3T positron emission tomography-MR scanner (Biograph mMR; Siemens Healthcare, Erlangen, Germany). The proposed approach was validated in a kiwifruit phantom, as its zones resemble the zonal anatomy of human prostate in T 2 -weighted (T 2 W) and DWI. 28 Feasibility for prostate ADC mapping was tested in three healthy subjects (mean age 33 ± 2 years). Local institutional review board approval and informed consent were obtained for all subjects.

| Acquisition sequence
We used a prototype 3D msDP sequence with gradient echo readout (msDP-GRE), shown in Figure 1. Each shot consisted of a collection of gradient echoes preceded by a twice-refocused spin echo module (T 2 -preparation [T 2 P]) with adiabatic pulses for improved robustness to B 0 and B 1 inhomogeneities. 29,30 DP was achieved by placing a pair of bipolar diffusion gradients symmetrically around the refocusing pulses to mitigate eddy currents. [22][23][24] Fat signal suppression was achieved using a water-selective 1-1 binomial excitation pulse. Data were acquired using a customdesigned 3D Cartesian SEctor-like Centric Trajectory with center OveRsampling (SECTOR; Figure 1). In each shot, the k-space points in k y -k z were sampled using a sector-like interleaf that can be undersampled in the periphery to speed up the acquisition (while 30% of the k-space center remains fully sampled). The key features of the SECTOR trajectory are: (1) centric acquisition, enabling efficient diffusion contrast encoding and minimization of T 1 effects; (2) k-space center oversampling at every shot; and (3) high number of shots, with features (2) and (3) enabling efficient magnitude averaging throughout the acquisition.
F I G U R E 1 Diagram of the proposed 3D multishot diffusion prepared-gradient echo sequence and SECTOR trajectory. The diffusion preparation (DP) module applied at the beginning of each shot is shown in more detail in the left box. The custom 3D SECTOR trajectory samples the k y -k z k-space plane on a Cartesian grid following sector-like interleaves. The first 10 segments of each GRE readout sample the same central k-space region (Kc) following a spiral ordering; the remaining GRE segments sample the k-space periphery (Kp) with the shot-specific sectorinterleaf. GRE, Gradient echo; RF, radiofrequency; SECTOR, SEctor-like Centric Trajectory with center OveRsampling; TR, pulse repetition time

| Acquisition parameters
Common msDP-GRE acquisition parameters for phantom and in vivo experiments were: 78 radiofrequency (RF) pulses per shot (10 of which were used for center oversampling), shot duration (TR) = 1000 ms, diffusion preparation-echo time (DP-TE) = 90 ms (to match the TE of the clinical EPI protocol), GRE FA = 10°, GRE TR/TE = 5.8/2.9 ms, 400 Hz/px bandwidth, transverse orientation, acquisition matrix 256 × 168 × 30. The resolution in the phantom was 1.0 × 1.0 × 5 mm 3 with field of view (FOV) = 260 × 170 × 150 mm 3 , and the resolution in vivo was 1.6 × 1.6 × 5 mm 3 with FOV = 420 × 278 × 150 mm 3 . Two b-values (50, 800 s/mm 2 ) and three diffusion directions (slice, phase, read) were acquired sequentially. Each volume was undersampled by a factor of two, resulting in an acquisition time per volume of 1 min (number of shots, N shots = 59) and a total acquisition time of 12 min. Three additional T 2 W volumes (T 2 P-TE = 0, 90, 150 ms) were acquired to generate a T 2 map needed for the ADC dictionary matching. 31 The T 2 -map acquisition was also undersampled by a factor of two, for an acquisition time of 3 min.

| Challenge 1: magnitude errors
To address the ME problem, a pair of phase-cycled volumes was acquired per b-value and diffusion direction, with the DP tip-up pulse of the second volume-phase shifted by 90°. For each phase-cycled volume, the high number of oversampled N shots k-space centers was averaged to stabilize the intershot MiME. Each averaged undersampled k-space was reconstructed using tvSENSE using MATLAB (MathWorks, Natick, MA), and the final corrected volume was calculated via square root sum of squares (SRSS) of the two phasecycled magnitude reconstructions to correct for EiME.

| Numerical simulations
A framework was developed to simulate the MEs in multishot acquisitions, as depicted in Figure 2. Spatially varying MEs can be simulated drawing random values from a normal distribution of phase errors, where the mean of the normal distribution corresponds to the shot-independent EiME phase F I G U R E 2 Framework to simulate the magnitude errors (MEs) in a multishot acquisition (simplified scenario with N shots = 5). A, For each shot i, a phase-error distribution is generated with phase offset and SD Δϕ MiME , representing eddy current-induced magnitude errors (EiMEs) and motion-induced magnitude errors (MiMEs), respectively. B, The correspondent magnitude error map, ME i is calculated as the cosine of the phase map (or the sine when simulating a 90° phase shift for phase cycling). C, The numerical phantom I reference is multiplied by the shot-specific ME i . D, The corrupted phantom is Fourier-transformed, and the k-space sampled with the correspondent shot i of the trajectory. E, The steps (A-D) are repeated N shots times, each time drawing different random values from the normal distribution to simulate intershot MiME variations. F, The corrupted N shots are then combined in a single multishot k-space averaging the oversampled k-space centers, finally obtaining the corrupted phantom I corrupted | 1445 ROCCIA et Al.

offset (
) and its SD to the intershot MiME-phase variations (Δϕ MiME ).
To simulate realistic MEs, magnitude errors observed in an in vivo prostate msDP-GRE acquisition were estimated (ME prostate ), and then mimicked using the simulation framework given in Figure 2. The ME prostate were estimated as follows: (1) for each shot i of the msDP-GRE acquisition (i = 1: N shots ), a low-resolution image was reconstructed from the oversampled k-space center; (2) a template was selected in a volumetric region of interest (ROI) around the prostate; and (3) the mean signal intensity of this template was calculated and plotted as a function of shots ( Figure 3). The MEs observed in a static kiwifruit phantom acquisition were also calculated for comparison with the in vivo ME prostate .
A numerical phantom (I reference ) was generated in MATLAB to mimic pelvic anatomy, 32 and was corrupted (I corruputed ) with the simulated ME as seen in Figure 2. A new set of magnitude errors (ME +90 ) was then simulated to generate the correspondent phase-cycled image (I corruputed+90 ). The ME and ME +90 pair was simulated using the same normal distribution, but drawing two different sets of random values to mimic a realistic phase-cycling scenario where the EiME is constant, but the intershot MiME of the twophase-cycled images may differ. The proposed k-space center averaging and phase-cycling approach were then used to recover the signal magnitude (I SRSS ). The signal intensity of I reference , I corrupted , I corruputed+90 , and I SRSS was compared in two ROIs: simulated transition zone (TZ) and peripheral zone (PZ). To simulate measurement noise in the numerical phantom, the experiment was repeated adding random white Gaussian noise in k-space, 33 resulting in an SNR = 14.5 and 11.1 for TZ and PZ, respectively. Statistical differences were assessed using a paired-sample Student t test with threshold P = .05.

| Phantom experiment
The proposed k-space center averaging and phase-cycling strategy was validated in the kiwifruit phantom with respect to EiMEs, as the MiMEs are negligible in a static phantom. Because a different amount of EiMEs is expected for each different diffusion direction, the signal-intensity dependence on slice, phase, and read directions was evaluated within a ROI to highlight the different amount of signal loss experienced in each direction. The signal dependence on diffusion direction in the corrupted and SRSS-corrected msDP-GRE was then compared with the ssDW-EPI signal diffusion dependence.

| Challenge 2: T 1 effects
To address challenge 2, the EPG formalism was used to model the acquisition-specific diffusion signal evolution along the readout and to generate a dictionary of simulated signals for ADC dictionary matching. Indeed, a simple monoexponential fitting would not take into account the consequences of these T 1 -related effects on the magnetization, and therefore would not accurately describe the msDP-GRE signal decay, as shown in Supporting Information Figure S1. A dictionary of signals was generated as a function of b-value, averaging the first 50% of the segments of each readout, which dominates the image contrast encoded with the SECTOR trajectory. To determine quantitative ADC values, F I G U R E 3 Simulated magnitude errors (MEs) (A) that mimic the errors measured in a prostate multishot diffusion prepared-gradient echo acquisition at b-value = 800 s/mm 2 (B). The dotted lines represent the intershot motion-induced variations (motion-induced magnitude errors [MiMEs]), and the solid lines represent the signals obtained averaging the intershot MiMEs (magenta for ME and yellow for ME +90 ). The deviation of the obtained signals from the reference one (solid-black line) represents the eddy currents magnitude error. The solid-green line is the square root sum of squares (SRSS) of the two phase-cycled averages, which is compared to the reference signal in the simulated scenario (A), with only a 0.6% signal difference. The MEs measured in the phantom, plotted in (C), are considerably smaller than the in vivo errors (B) as the phantom is static after normalizing both acquired data and simulated dictionary by the respective first data point (corresponding to b-value = 50 s/mm 2 ), the acquired diffusion data were matched via voxel-wise L 2 -norm minimization to the dictionary.

| Numerical simulations
To characterize the dependence of the diffusion signal decay on T 1 and T 2 relaxation times, a set of signals extracted from the simulated dictionary was plotted for a range of T 1 = [1600:100:2000] ms and T 2 = [80:20:160] ms. Simulations were performed for typically reported values in healthy and cancerous prostate PZ at 3T (PZ healthy : T 2 /ADC = 160 ms/ 1.6 × 10 −3 mm 2 /s, PZ cancer : T 2 /ADC = 100 ms/1.0 × 10 −3 mm 2 /s). 15,[34][35][36][37] The sensitivity of the signal to T 1 and T 2 and the PZ healthy versus PZ cancer contrast was calculated as the percentage difference of the respective signals.

| Phantom and in vivo experiments
Although the proposed simulation-based approach takes into account T 1 effects caused by TR < T 1 and by T 1 regrowth during readout, it does not take into account the T 1 effects caused by DP tip-down/up FA imperfections. Therefore, a phantom experiment was performed repeating the msDP-GRE acquisition with different effective FAs, within a range of expected variations: FA = [80, 85, 88, 90, 92, 95, 100]° ( Figure 6C).
The proposed dictionary-matching approach was then validated in the kiwifruit phantom. The acquired images at b-values = 50 and 800 s/mm 2 were first reconstructed using the ME-correction method described above. ADC maps were then obtained with both conventional monoexponential fitting and dictionary matching. Based on the simulation results, the dictionary matching used a fixed T 1 value of 1700 ms (representative of prostate T 1 at 3T) and a voxel-specific T 2 map. To evaluate the impact of this choice on the ADC estimation, a sensitivity experiment was performed assuming different fixed T 1 and T 2 values in the ADC matching (Supporting Information Figure S2). The estimated ADC values were compared with clinical standard ssDW-EPI ADC values in the outer pericarp (OP) and peripheral placenta (PP) of the kiwifruit, as these zones have properties similar to prostate PZ and TZ. 28 Statistical difference was assessed using a paired-sample Student t test with threshold P = .05 (Figure 7).
Feasibility of the ADC dictionary matching was finally tested in three healthy subjects. Images acquired with the proposed msDP-GRE sequence were corrected for MEs and reconstructed with tvSENSE. Geometric distortion was evaluated by manually drawing contours around the prostate gland and rectum. The Dice similarity score 38 was then calculated between the ssDW-EPI and msDP-GRE contours and those obtained on the reference anatomical T 2 W-TSE. The sensitivity of the ADC matching to T 1 and T 2 variations was also assessed for the prostate scan (Supporting Information Figure S2).

| Numerical simulations
The simulated ME and phase-cycled ME (ME +90 ) are plotted as a function of shot in Figure 3A. These MEs were generated to mimic the MEs estimated in vivo (ME prostate ; Figure 3B), with an even slightly higher variation (24% and 17% for ME and ME prostate , respectively). Such ME distribution was achieved simulating a phase-error distribution with offset = 46° and SD Δϕ MiME = 10°. The average of ME and ME +90 describes the effect of averaging the oversampled k-space center throughout the shots. The SRSS correction applied to the average of the simulated ME and ME +90 successfully recovered the magnitude loss; indeed, the SRSS signal matched the reference uncorrupted signal magnitude with only a 0.6% magnitude difference ( Figure 3A). Figure 3C shows the ME estimated in the phantom, where the intershot MiMEs were considerably lower than those observed in vivo as the phantom is static. The EiMEs (offset in MEs and ME +90 ) were still present in the phantom as these depend on the system. The magnitude of the numerical abdominal phantom corrupted with the MEs of Figure 3A could be recovered using the proposed k-space center averaging and phase-cycling approach ( Figure 4). The signal intensity measured in the corrected phantom (I SRSS ) for two ROIs (TZ and PZ) showed no statistically significant difference compared to the reference phantom (I reference ). When adding noise (SNR = 14.5 and 11.1 for TZ and PZ, respectively), the signal intensity could still be recovered, but the SD was higher: from 0.042 to 0.057 for TZ, and from 0.022 to 0.053 for PZ.

| Phantom experiment
The acquired msDP-GRE phase-cycled images in the phantom are shown in Figure 5, together with the SRSS correction (msDP-GRE SRSS ), for both b-values = 50 s/mm 2 and 800 s/mm 2 . For the low b-value (50 s/mm 2 ), the msDP-GRE manifested little magnitude corruption for all diffusion directions, indicating that the EiMEs were very mild; indeed, the 90° phase-cycled images (msDP-GRE +90 ) contain almost no signal. This is confirmed with the ROI analysis depicted in Figure 5C, which shows also that the (minor) SRSS | 1447 ROCCIA et Al.
correction led to a signal intensity with similar behavior to ssDW-EPI across all diffusion directions. The SRSS correction was able to recover the magnitude for all directions, and the dependence of msDP-GRE SRSS on the diffusionsensitizing gradient direction observed within the ROI was also observed with the reference ssDW-EPI ( Figure 5D).

F I G U R E 4
Reference numerical phantom (I reference ), corrupted phase-cycled images I corrupted , and I corrupted+90 (obtained applying magnitude error (ME) and ME +90 of Figure 3A), and phantom corrected via square root sum of squares (SRSS) of I corrupted and I corrupted+90 (I SRSS ). No statistically significant difference in signal intensity was observed between I reference and I SRSS , and also between I reference and I SRSS + noise , where Gaussian noise was added to the image prior to magnitude corruption (signal-to-noise ratio = 14.5 and 11.1 for transition zone and peripheral zone, respectively). ns, Nonsignificant. *P < .05 F I G U R E 5 Kiwifruit phantom multishot diffusion prepared-gradient echo (msDP-GRE) acquisition, with the two phase-cycled images (+90° indicating the phase shift) and the square root sum of squares (SRSS)-corrected images at b-values = 50 s/mm 2 (A) and 800 s/mm 2 (B). The eddy-current-induced magnitude corruption can be observed especially at the high b-value in the slice and read directions, where a phase offset causes the signal to lie in both phase-cycled images (so on both axes of the transverse plane). The bar plots in (C) and (D) show the signal intensity dependence on the diffusion direction within a region of interest, and is compared to the reference single-shot diffusion-weighted echo planar imaging (ssDW-EPI) direction dependence. The direction-dependent variations observed with the corrected msDP-GRE SRSS are comparable to those observed with ssDW-EPI

| Numerical simulations
The simulated signal evolution for the proposed sequence as a function of T 1 and T 2 is shown in Figure 6A,B. The percentage signal difference between PZ healthy and PZ cancer was 36%, underlying desired ADC contrast. Only slight variations were observed over the range of simulated T 1 values ( Figure 6A), with a 3% signal difference between the T 1 = 1600 ms and T 1 = 2000 ms signal, suggesting that the T 1 sensitivity is negligible compared with the observed diffusion contrast (36%). This provides a rationale for using a global fixed T 1 (T 1 = 1700 ms) for all voxels in the dictionary matching. On the other hand, the diffusion signal is sensitive to T 2 variations ( Figure 6B), with a 17% and 21% signal difference between the simulated T 2 = 80-ms signal and T 2 = 160-ms signal for PZ healthy and PZ cancer , respectively. This finding suggests that different T 2 species require a different dictionary entry to describe their signal behavior, so a T 2 map (calculated from the three T 2 W volumes via dictionary matching 31 ), should be included in the matching to guide the dictionary search and estimate the correct ADC for each tissue type.

| Phantom and in vivo experiments
Signal variations measured in the kiwifruit phantom ROI caused by imperfect DP tip-down/up pulse FAs are shown in Figure 6C. The percentage difference of the measured signal across the range of effective FAs was 3.6% and 6.9% for bvalue 50 and 800 s/mm 2 , respectively, demonstrating robustness of msDP-GRE to FA variations.
For the healthy subjects' msDP-GRE acquisition, the ME correction and tvSENSE reconstruction provided F I G U R E 7 Kiwifruit apparent diffusion coefficient (ADC) maps obtained with single-shot diffusion-weighted echo planar imaging (ssDW-EPI) and monoexponential fitting, and proposed multishot diffusion prepared-gradient echo (msDP-GRE) using both conventional monoexponential fitting and dictionary matching. The estimated ADC is compared in the kiwifruit outer pericarp (OP, shaded green area) and peripheral placenta (PP, shaded yellow area), as these zones have similar tissue properties to prostate peripheral zone and transition zone. 28 No statistical difference was observed between ADC values obtained with clinical standard ssDW-EPI and proposed msDP-GRE with dictionary matching, whereas ADC values obtained using msDP-GRE with monoexponential fitting were significantly different from ssDW-EPI ADC. The ADC values could not be accurately estimated at the center of the kiwifruit because of the short T 2 in this region (T 2 = 45 ms), which is lower than T 2 values reported for prostate cancer (T 2 > 70 ms [35][36][37] ); therefore, this zone was masked out. ns, Nonsignificant. *P < .05 F I G U R E 8 Geometric distortion analysis for a healthy-subject prostate acquisition (b-value = 50 s/mm 2 ). Outline of prostate and rectum are drawn (in solid and dotted lines, respectively) on reference anatomical T2-weighted-turbo spin echo (T 2 W-TSE; A, yellow line), clinical standard single-shot diffusion-weighted echo planar imaging (ssDW-EPI; B, magenta line), and proposed multishot diffusion prepared-gradient echo (msDP-GRE; C, blue line). Outlines from A, B, and C are overlapped in D to compare geometric distortion. The msDP-GRE outlines matched the reference T 2 W-TSE, whereas the ssDW-EPI outlines deviated more from the reference good-quality images. The msDP-GRE images showed no geometric distortion when compared with the T 2 W-TSE, whereas in the ssDW-EPI images moderate geometric distortion was observed at the air-tissue interface, as highlighted by the anatomical outline comparison in Figure 8. The proposed msDP-GRE method achieved a higher mean Dice similarity score compared with ssDW-EPI in the prostate (0.91 ± 0.02 and 0.86 ± 0.02, respectively) and rectum (0.74 ± 0.01 and 0.64 ± 0.07, respectively). The Dice similarity scores for all healthy subjects are listed in Table 1.
A representative case of an ADC map obtained with the proposed sequence and dictionary matching is shown in Figure 9, alongside the clinical standard ssDW-EPI ADC map and the dictionary-based T 2 map used to guide the ADC dictionary matching. The estimated ADC in the PZ was 1.3 ± 0.18 × 10 −3 mm 2 /s, in agreement with reported literature values in young healthy subjects. 39,40

| DISCUSSION
In this study, we proposed a framework for distortion-free 3D prostate ADC mapping that combines a msDP acquisition with ME compensation, tvSENSE reconstruction, and ADC-dictionary matching. A simulation framework was developed to describe and mimic the ME corruption problem in multishot acquisitions. The proposed approach was assessed in simulations and phantom experiments, and initial in vivo feasibility was demonstrated in healthy subjects.
Obtaining distortion-free diffusion images is crucial in prostate imaging. Susceptibility-induced distortions at the rectal interface may hamper the use of diffusion imaging in applications requiring high 3D geometric accuracy, such as high-precision diagnostic and MR-guided radiotherapy planning. Here, with the proposed 3D msDP-GRE sequence, we achieved better geometric fidelity than ssDW-EPI using a GRE readout.
The multishot sequence holds potential to achieve high spatial resolution that, together with the 3D acquisition, may provide a comprehensive characterization of the whole prostate, and an accurate detection of small cancerous lesions that could otherwise be missed. 41 In msDP sequences, the two main challenges that cause inaccurate ADC quantification are MEs caused by intershot motion and eddy currents (MiMEs and EiMEs), and T 1 effects, both confounding the desired diffusion weighting.

| Challenge 1: magnitude errors
The first challenge was addressed using bipolar gradients, [22][23][24] and a combination of oversampled k-space center averaging and phase cycling. Other studies have previously used phase cycling to account for motion and eddy currents in DP sequences, 17,25 at the cost of doubling the scan time.
Here, each phase-cycled acquisition was undersampled by a factor of two, ie, the total acquisition time was not increased.
The oversampling property of the SECTOR-trajectoryenabled acquisition and averaging of the k-space center at T A B L E 1 Dice similarity score for the clinical standard ssDW-EPI and the proposed msDP-GRE (b-value = 800 s/mm 2 ) compared with the ground truth T 2 W-TSE, for prostate gland and rectum of all healthy subjects F I G U R E 9 Representative case of a healthy subject prostate acquisition: T 2 map used in the apparent diffusion coefficient (ADC) dictionary matching (A), proposed nondistorted three-dimensional ADC map obtained with multishot diffusion prepared-gradient echo (msDP-GRE) and dictionary-based matching (B), and clinical standard two-dimensional single-shot diffusion-weighted echo planar imaging (ssDW-EPI) ADC map (C) every shot. This, together with the high number of shots, helped recovering intershot magnitude instabilities and thus magnitude corruption when combining the two-phase-cycled images with SRSS. This was validated first in a numerical phantom experiment, where realistic intershot MEs were deliberately introduced using the proposed simulation framework, and the signal magnitude could be recovered using the proposed approach. In the kiwifruit phantom experiment, the EiME effect was much stronger when increasing the b-value, particularly in the slice and read directions where the signal was mostly sampled by the extra phase-cycled acquisition, indicating the presence of an eddy current-induced-phase offset. The remaining intershot magnitude variations measured in the static phantom ( Figure 3C) can be explained with table vibrations caused by fast-switching diffusion gradients. The SRSS reconstruction enabled the recovery of the signal magnitude for all diffusion directions. The k-space oversampling property of the proposed SECTOR trajectory for increased robustness to intershot motion is similar to PROPELLER (periodically rotated overlapping parallel lines with enhanced reconstruction) 42 or variable-density spirals 43 methods. However, the proposed method has the main advantage of acquiring the data in 3D, and the SECTOR trajectory samples the k-space points on a Cartesian grid so it does not require sampling density correction and gridding steps during reconstruction. In addition, combining this type of msDP-GRE acquisition with a dictionary-matching approach enables taking into account potential modification of the trajectory, eg, the extent of oversampled center by simply averaging a different amount of readout segments when generating the dictionary.

| Challenge 2: T 1 effects
The influence of the T 1 regrowth was mitigated in the first place using a centric trajectory, which immediately encodes the DP magnetization in the center of k-space. 18,19,24 Then, the T 1 -related effects (incomplete T 1 recovery and T 1 regrowth) that make the msDP-GRE decay deviate from a monoexponential decay, were taken into account through the dictionary-based ADC matching, which simulates the acquisition-specific timing parameters and magnetization evolution, thus accurately describing the msDP-GRE signal decay.
The simulation results have shown that T 1 variations as typically found in the prostate, in combination with the centric trajectory, provokes only small variations in the msDP-GRE diffusion-decay signal (Figure 6A), providing a rationale for using a global T 1 value for all voxels in the ADC matching, similar to the findings in our previous study. 31 On the other hand, the ADC matching is sensitive to T 2 variations, so the underlying voxel-specific T 2 values need to be included in the ADC-matching routine. This was also demonstrated in phantom and prostate scans (Supplementary Information Figure S2).
The phantom experiment demonstrated robustness of the proposed msDP-GRE sequence to effective FA variations of the DP tip-down/tip-up pulse. The ADC values in the phantom experiment could be accurately estimated using msDP-GRE and dictionary matching, whereas significant ADC bias was observed when using the monoexponential fitting, confirming its inadequacy to accurately describe all the effects occurring in the magnetization evolution. A higher SD was observed in msDP-GRE ADC maps compared with ssDW-EPI, likely because of ssDW-EPI blurring and partial-volume effects, interpolation, and smoothing. The geometric distortion analysis demonstrated good performance of the proposed msDP-GRE sequence. The Dice similarity score for the rectum was lower than for the prostate because of rectal distension and contraction between the scans. The ADC measured in the prostate PZ was in agreement with literature values in young healthy subjects. 39,40

| Study limitations
It should be noted that intershot magnitude variations caused by physiological motion expected in the prostate are rather benign, as shown in Figure 3B, and as already discussed in other studies. 25 The robustness of the proposed approach in organs with substantial physiological motion (eg, brain, heart) remains to be investigated.
Because the dictionary matching uses a fixed global T 1 value for all the voxels, it does not account for structures with T 1 values outside the range of negligible T 1 sensitivity found in this study. To address this limitation, a T 1 map, which is typically acquired in the clinical mpMRI prostate protocol for DCE analysis, could be included in the dictionary matching.
Another limitation of the current framework is that the ADC could not be correctly estimated in the kiwifruit center, which is characterized by short T 2 -relaxation times and low SNR. Indeed, with short T 2 -relaxation times, the DP signal decays very quickly, while the SECTOR trajectory is still sampling the k-space center, so that the signal is dominated by T 1 contrast from magnetization recovery towards the GRE steady state, prohibiting an accurate diffusion-contrast encoding. Although this effect could generally be accounted for using dictionary matching, this was not possible for such low T 2 (mean T 2 = 45 ms). The T 2 of the kiwifruit center was lower than reported T 2 -relaxation times for PZ cancerous tissue. [35][36][37] Optimization of the SECTOR trajectory will be investigated in future work to achieve a trade-off between number of segments used to oversample the center for intershot motion averaging, and diffusion contrast encoding. An optimized variable FA scheme will also be investigated in future work to address the T 1 regrowth issue.
In this work, the healthy subject population was small, and clinical validation in patients still needs to be performed. For clinical acceptance, the scan time efficiency of the proposed sequence should be improved, eg, by further increasing the acceleration factor, exploiting the flexible undersampling capability of the trajectory in the two phase-encoding directions.

| CONCLUSION
The proposed msDP-GRE sequence together with dictionary matching produced 3D distortion-free ADC maps of the prostate. The oversampled k-space center averaging and phase-cycling strategy was able to recover signal magnitude loss, and the simulation-based ADC dictionary matching was used to account for T 1 effects. The proposed approach may be advantageous when high 3D geometric accuracy is desired. Future work will entail optimization of the trajectory, and a larger study on healthy subjects and prostate cancer patients for clinical validation.