Repeatability of IVIM biomarkers from diffusion-weighted MRI in head and neck: Bayesian probability versus neural network.

PURPOSE: The intravoxel incoherent motion (IVIM) model for DWI might provide useful biomarkers for disease management in head and neck cancer. This study compared the repeatability of three IVIM fitting methods to the conventional nonlinear least-squares regression: Bayesian probability estimation, a recently introduced neural network approach, IVIM-NET, and a version of the neural network modified to increase consistency, IVIM-NETmod .
METHODS: Ten healthy volunteers underwent two imaging sessions of the neck, two weeks apart, with two DWI acquisitions per session. Model parameters (ADC, diffusion coefficient D t , perfusion fraction f p , and pseudo-diffusion coefficient D p ) from each fit method were determined in the tonsils and in the pterygoid muscles. Within-subject coefficients of variation (wCV) were calculated to assess repeatability. Training of the neural network was repeated 100 times with random initialization to investigate consistency, quantified by the coefficient of variance.
RESULTS: The Bayesian and neural network approaches outperformed nonlinear regression in terms of wCV. Intersession wCV of D t in the tonsils was 23.4% for nonlinear regression, 9.7% for Bayesian estimation, 9.4% for IVIM-NET, and 11.2% for IVIM-NETmod . However, results from repeated training of the neural network on the same data set showed differences in parameter estimates: The coefficient of variances over the 100 repetitions for IVIM-NET were 15% for both D t and f p , and 94% for D p ; for IVIM-NETmod , these values improved to 5%, 9%, and 62%, respectively.
CONCLUSION: Repeatabilities from the Bayesian and neural network approaches are superior to that of nonlinear regression for estimating IVIM parameters in the head and neck.

INTRODUCTION

Magnetic resonance DWI is used for diagnostic and prognostic purposes in head and neck cancer. , , , In DWI, signal decreases with diffusion weighting as result of Brownian motion of water molecules and other intravoxel incoherent motions (IVIMs) (ie, “microscopic translational motions that occur in each image voxel”). , By fitting the DWI signal from different diffusion weightings to an exponential model, its parameters can be estimated. A mono‐exponential can be used to estimate the ADC. A bi‐exponential model (the IVIM model ) can be used to additionally model pseudo‐diffusion component () and perfusion fraction ()—both related to the microcirculation of blood—resulting in the corrected or “true” diffusion coefficient (). Because the restriction of diffusion is related to the microstructure of tissue (eg, cellular density), this can characterize tumors and provide early information on changes due to (or despite) treatment occurring before detectable tumor growth or shrinkage.

The IVIM model is appealing, as it allows the assessment of the additional biomarkers and . However, IVIM parameter estimation tends to be very sensitive to noise. As a result, parametric maps are often noisy and show poor repeatability. Poor repeatability limits the use of IVIM in practice, because precision is required for patient‐specific clinical use of IVIM.

Recently, novel fitting methods with a Bayesian probability approach , , and a neural network have shown promising results in terms of reduced noise in the parameter maps based on simulations, and they reduced interobserver variability in vivo. If these techniques also help improve test–retest repeatability in vivo, they could help introduce IVIM into clinical workflows.

Therefore, in this study we investigate these new methods in terms of test–retest repeatability. We compare the intrasession and intersession repeatability of the least‐squares fitting method, the Bayesian inference fitting method, and two neural network–based fitting methods for in vivo IVIM data in the head and neck region in healthy volunteers. We hypothesize that the new Bayesian and neural network approaches will outperform the conventional least‐squares fitting approach.

METHODS

This study was approved by the local medical ethics committee, and written informed consent was obtained from all subjects. Ten healthy volunteers were included: 7 males, 3 females, mean age 33 years (range 22‐50 years). Each volunteer underwent two MRI sessions (at least 2 weeks apart) with two examinations per session. The subject was taken out of the MR scanner between examinations. Sequences were acquired on a 3T Ingenuity TF PET/MR scanner (Philips Healthcare, Best, the Netherlands) equipped with a 16‐channel neurovascular coil. Each examination consisted of an axial stack of 29 T1‐weighted turbo spin‐echo images followed by a stack of DWI acquisitions in the same 29 imaging planes, covering the neck from the larynx until the base of the skull. Diffusion‐weighted imaging was acquired with a single‐shot spin‐echo EPI sequence with 12 b‐values (0, 2, 5, 25, 50, 75, 100, 150, 300, 500, 700, and 1000 s/mm2). Only the DWI images with b = 1000 s/mm2 were averaged over two acquisitions. Diffusion weighting was performed in three orthogonal directions with bipolar gradients, TE = 57 ms, TR = 3242 ms, gradient time interval = 28 ms, and gradient duration = 18 ms. Further scan parameters were as follows: acquired matrix size = 128 × 111 × 29, acquired voxel size = 1.88 × 1.95 × 4 mm3, reconstructed voxel size = 1 × 1 × 4 mm3, and short TI inversion recovery used for fat suppression with a 230‐ms TI. The DWI scan duration was 6 minutes. Motion correction of the DWI images was applied by image registration, as provided by the scanner software.

The DWI data were processed voxelwise to generate parametric maps of the ADC and IVIM parameters. Parameter estimates were extracted for two tissues: tonsil and medial pterygoid muscle. Volumes of interest (VOIs) were defined on the images without diffusion weighting ( s/mm2) while using the T1‐weighted image for anatomical reference. The T1 images were not co‐registered to the diffusion‐weighted images. Delineation was performed using in‐house‐developed software by a single observer in one session. Spherical VOIs of 5‐mm radius were placed in each tonsil, and spherical VOIs of 6‐mm radius were placed in the medial pterygoid muscle on each side. These VOIs were small enough to always fit inside the tissues of interest. The VOIs were projected onto the parametric maps, all voxels with (partial) overlap were extracted, and the median values of the parameters were then calculated.

The signal at s/mm2 was excluded (except for calculating fit boundary of [see subsequently] and for normalization purposes in the neural network) for the parameter estimations described subsequently, reducing the number of b‐values to 11. The reason for this is to reduce attenuation effects of macroscopic flow at small b‐values. , This additional accelerated decay between s/mm2 and the first nonzero b‐value is not accounted for in the conventional IVIM and ADC models. ,

The mono‐exponential model used to estimate the ADC is given bywhere is the signal intensity without diffusion weighting ( s/mm2), and and are estimated by performing a linear least‐squares fit on the log‐transformed data, as implemented by the scanner manufacturer.

The IVIM model extends the ADC model with a second exponential. The bi‐exponential equation of the model is given bywhere is the perfusion fraction; the pseudo‐diffusion coefficient; is the diffusion coefficient; and is the fitted signal intensity for  s/mm2. The IVIM model parameters were estimated using four different approaches: a nonlinear least‐squares fit, a Bayesian approach, and two neural network‐based fitting approaches. The two neural network approaches consisted of a network nearly identical to the original publication (IVIM‐NET) and a modified network (IVIM‐NETmod), as detailed later.

Nonlinear least squares

The nonlinear least‐squares fit was performed using the trust‐region reflective algorithm as implemented in MATLAB R2019a (MathWorks, Natick, MA), with the following fit boundaries: , mm2/s, mm2/s, and . Starting values were selected randomly in the range within the fit boundaries as provided by MATLAB functionality.

Bayesian probability

The Bayesian approach was also performed in MATLAB R2019a and was based on a previous publication. In short, the method gives a maximum a posteriori estimate of each parameter by maximizing the marginal posterior probability density functions, which are acquired by means of slice‐sampling the joint posterior probability. , , The following multiparametric Gaussian likelihood function was used:where is the number of b‐values. The constraint was implemented in the joint prior distribution. Lognormal distribution priors were used for and ; a beta distribution prior was used for ; and a uniform distribution prior was used for . The priors for , , and were estimated by fitting these distributions to results of a prerun of the same Bayesian approach using bounded uniform priors (, mm2/s, mm2/s, and ).

Neural network

The IVIM‐NET approach was carried out in Python 3.7.4 and PyTorch 1.3.0 using the open access code from the original publication , and code obtained from the repository of co‐author OGC (currently shared on request; will be made public in the near future and maintained as IVIM‐NET evolves). Our source code with the network definitions and training methods is available on GitHub.

The network, depicted in Figure 1, consists of an input layer with a number of neurons equal to the number of b‐values used to analyze the data, three fully connected hidden layers (each with the same number of neurons, each using the exponential linear unit activation function), and an output layer with a neuron for each parameter. Background voxels were excluded by manually thresholding the  mm2/s images. Training was performed on the entire data set for each epoch, combining and shuffling the voxels from all patients. Data normalization, which is standard for neural networks, was performed using S(0). This signal measured at s/mm2 was only used for normalization, not as input for the network. The mean squared error between the fitted and actual, normalized, signal was used as loss function. An early stopping criterium (patience) of 10 bad epochs was used, so training was stopped when no improvement was found during the last 10 epochs. Different from the original publication, we included an output neuron for , where and are the estimated and measured signal intensity at  s/mm2, respectively.

FIGURE 1

Neural network architecture, created with NN‐SVG. The network predicts x1 to x4, which are converted to the intravoxel incoherent motion (IVIM) parameters by the constrain function g(x) using Equations 4 (original network) and 5 (modified network), to add parameter constraints

In addition to this implementation, we made a few modifications in a new implementation, IVIM‐NETmod. The IVIM parameters were constrained by g(x). In the original network, the predicted IVIM parameters were constraint by taking the following absolute:

In the presented modified network, a sigmoid function was applied to the output as constraint instead:which rescaled the output between the following fit boundaries (min < parameter < max): , mm2/s, mm2/s, and . Second, with the aim of preventing overfitting, we split the data set into two parts: one for training (80%) and one for validation (20%). For the same reason, we reduced the patience (early stopping criterion, see previously) from 10 to 4. Furthermore, as we had a substantially larger data set than Barbieri et al, we limited the number of iterations during each training epoch to 1024, such that we regularly validate how well the network is performing even for large data sets. Because the batch size of an iteration is fixed (128 voxels), each epoch no longer processes the entire data set but a random selection of the training set (in our case, approximately 1.5%). Each epoch evaluates the entire validation set.

Network consistency

To investigate the consistency of the IVIM‐NET approaches as a whole (ie, whether the network converges to consistent estimates), we repeated the complete process of training 100 times. Each time, the network was initialized with new random weights and shuffling (and splitting in the case of IVIM‐NETmod) of the data set. We compared the runs qualitatively by visual inspection of the parametric maps. We investigated consistency by calculating the average parameter values for both tissues over all subjects and sessions in each run, and then calculated the coefficient of variance (CoV) over the 100 runs.

Statistics

The intrasession repeatability was calculated by considering the two measurements within a session as paired measurements. Conversely, the intersession repeatability was calculated by considering the first measurements in each session as one pair, and the second measurements in each session as another pair. Moreover, the left and right measured values were considered measurements for the same tissue of interest (ie, tonsil and pterygoid muscle). Thus, each subject had four pairs of observations for the calculation of repeatability, and pairs were either between sessions or within sessions.

We used 95% confidence intervals of the mean difference between paired measurements over all subjects (for both intrasession and intersession pairing) to verify that the repeated measurements were not systematically different. We then calculated the within‐subject coefficient of variation (wCV), which is a relative measure of repeatability. An overview of the concepts repeatability and consistency can be found in Table 1.

TABLE 1

Explanation of analysis concepts used in this study

Concept	Description	Quantification	Application
Repeatability	Variation between repeated measurements	wCV	All methods
Consistency	Variation between training runs of IVIM‐NET on same measurements	CoV	IVIM‐NET

We compared the repeatability of the four methods with paired Wilcoxon signed‐rank tests of the wCV estimates. For the IVIM‐NET methods, we calculated the median wCV of the 100 runs for each subject, and used these wCV estimates in the paired tests among the four methods. A P‐value below .05 was considered significant.

RESULTS

Figure 2 shows examples of the parametric maps calculated with the different methods. Parametric maps calculated by nonlinear regression were most noisy, followed by the Bayesian probability approach. The IVIM‐NET method showed the least noise and most anatomical detail, and was in these terms comparable to the ADC map. The maps estimated with nonlinear regression showed systematically higher values than the other methods, as shown in Figure 2. The maps of nonlinear regression and IVIM‐NETmod showed many regions with very high values.

FIGURE 2

Typical parametric maps of the estimated S0, ADC, true diffusion coefficient D. pseudo diffusion coefficient D, and perfusion fraction f. Regions of interest delineating the tonsils are shown in the ADC map. The pterygoids are not situated at this level; examples of regions of interest can be found in the Supporting Information. Parametric maps of other IVIM network (IVIM‐NET) instances can also be found in the Supporting Information. Abbreviation: IVIM‐NETmod, modified IVIM‐NET

None of the methods showed a systematic difference between the repeated measurements for any of the parameters. The calculated wCV estimates are shown in Figure 3. Notably, intrasession and intersession wCV was comparable and both VOIs show the same patterns when comparing the IVIM methods. The methods differ in terms of repeatability and, in general, wCV was highest (worst) when parameters were estimated using nonlinear least squares (except for in the pterygoid).

FIGURE 3

Within‐subject coefficient of variation (wCV) of the parameters for each method. The median value of 100 training runs is displayed for the neural network methods (*P ≤ 0.05, **P ≤ 0.01). Abbreviations: LLS, linear least squares; NLS, nonlinear least squares

This difference was often significant, especially for . Significance is indicated in Figure 3, comparing the wCV values of each of the methods for each of the parameters. Tables of the p‐values are available in Supporting Information Tables S1‐S3. The median repeatability results of IVIM‐NET and IVIM‐NETmod were mostly comparable to the repeatability of the Bayesian approach, except for , in which the wCV of IVIM‐NET was significantly better; IVIM‐NETmod was only significantly better for in the tonsils. The median repeatability of IVIM‐NET was better than for IVIM‐NETmod, although the difference was rarely significant.

Network consistency

Visual comparison of the parametric maps of the repeated network instances showed inconsistencies for both IVIM‐NET and IVIM‐NETmod; examples of this can be found in Supporting Information Figures S2‐S4. One hundred figures of the maps of each method are included in the Supporting Information. The network instances generally produced maps with a similar distribution but different offset/scaling values. Comparing the average parameter values for IVIM‐NET, the instances showed a CoV of 15% for and , and 94% for . The CoV for the pterygoids and tonsils were equal. The maps sometimes showed a visually different distribution. The wCV values for the IVIM‐NET instances were also inconsistent, as shown by the box plots in Figure 4. For IVIM‐NETmod, the average parameter values were more consistent, with CoVs of 5%, 9%, and 62% for , , and , respectively (again equal for pterygoids and tonsils). The wCV values for and were also more stable for IVIM‐NETmod, as reflected by the smaller intervals of the box plots in Figure 4.

FIGURE 4

Box plots (with Tukey‐type whiskers) of wCV values for 100 runs of the neural networks

DISCUSSION

In this study, we quantified the test–retest repeatability of nonlinear regression, neural network–based, and Bayesian IVIM in the head and neck region. Our results show that these latter two fit approaches substantially outperform the conventional nonlinear regression approaches commonly used for IVIM fitting. Furthermore, although IVIM‐NET has an improved test–retest repeatability, it has an additional uncertainty in that repeated training of networks gives inconsistent results on identical data.

Repeatability estimates of ADC in the tonsils using linear regression fit reported in this study are similar to those reported by Kang et al (also in the tonsils). Two other studies focused primarily on the primary lymph nodes, which makes it hard to compare results directly. Hoang et al report a repeatability coefficient in percentages (15%), which is equivalent to a wCV of 5.3%. The wCV reported by Paudyal et al is 2.38% and much lower than repeatability values found in this study. The generally larger volumes of metastatic lymph nodes might partly explain why the reported estimates are lower. In the case of the study of Paudyal et al, no repositioning of the subject in the magnet appears to have occurred between scans, which could be a major source of measurement variability. This might also explain the difference in reported wCV between the two studies. In our present study, the subject was taken out of the scanner between scans for the intrasession repeatability estimates. No differences were seen between intrasession and intersession repeatability. This indicates that long‐term (order of weeks) physiological variability over time was secondary to the measurement error and short‐term (~30 minutes) physiological variability.

The neural network approach for calculating IVIM maps was introduced only recently. Visual interpretation of the images suggests that more realistic parametric maps are produced by both neural networks compared with the other methods; these maps do not show isolated high or low pixels, and therefore appear to be least affected by noise in the acquisitions. This is in line with earlier observations from Barbieri et al. Our study has quantified the test–retest repeatability and shows that the network also outperforms linear regression regarding this aspect.

Network training for the entire data set took up to 1 hour for IVIM‐NET and up to 5 minutes for IVIM‐NETmod. Application of the network took only a couple of seconds for the entire data set. Barbieri et al had substantially less training data, and hence had training times of 5 minutes using the unmodified network. The large difference in training time in our data was primarily the result of decreasing the amount of data seen each epoch in the IVIM‐NETmod. This major advantage of analysis speed, compared with the other methods investigated in this study (around half an hour per scan with nonlinear regression, and multiple hours per scan using Bayesian probability fitting), makes it viable for use in clinical practice.

Although IVIM‐NET showed promising test–retest repeatability, consistency of the neural network approach is currently still an issue. Our results show that, after renewed training, the parameter values and repeatability estimates vary. The IVIM‐NETmod method showed more consistent results, although the method is still unstable for . Consistency of the approach might be improved by optimizing the starting point of the network, such as by choosing different weight initialization or by training on a set of simulated data first. Avoiding to fit (ie, fixing it instead to an a priori estimate) has been shown to improve repeatability , and might also improve network consistency.

A challenge for further research is to identify an acceptable neural network that does not only give estimates with good repeatability, but is also consistent after retraining. Until such consistency is achieved, it is imperative that a single network instance is used for comparative applications, such as in longitudinal studies. Use of separately trained networks will otherwise lead to biased results.

Although other DWI models , , are available, this study has been limited to the ADC and the IVIM model. Another limitation of our study is that we could not compare the methods in terms of accuracy, because a ground truth was unavailable in our study. We hope, therefore, that these methods will be included in future phantom studies. Finally, the choice of b‐values was probably not optimal; b‐value optimization may improve IVIM estimates. ,

CONCLUSIONS

The processing speed of the neural network makes it viable for use in clinical practice. However, the inconsistency of training results is challenging. Our presented modifications in the neural network make this approach more consistent, although the output still shows some inconsistency between different training runs on the same data set. Thus, the neural network approach needs to be further improved to identify neural networks that are both consistent and precise. Nonetheless, repeatability from the Bayesian and neural network approaches are superior to that of nonlinear regression for estimating IVIM model parameters.

FIGURE S1 The ADC maps with examples of the regions of interest indicated

FIGURE S2 Example of D maps from three instances of intravoxel incoherent motion network (IVIM‐NET) with similar distribution but different absolute values

FIGURE S3 Example of D maps from four instances of IVIM‐NET: two with similar distribution but different values (instances 1 and 2) and two with a different distribution (instances 5 and 11)

FIGURE S4 Example of D maps from three instances of a modified IVIM‐NET (IVIM‐NETmod) with different values

TABLE S1 Wilcoxon signed‐rank test p‐values, comparing paired within‐subject coefficient of variation (wCoV) values of the methods for D

TABLE S2 Wilcoxon signed‐rank test p‐values, comparing paired wCoV values of the methods for D

TABLE S3 Wilcoxon signed‐rank test p‐values, comparing paired wCoV values of the methods for f

Click here for additional data file.