Spectral diffusion analysis of kidney intravoxel incoherent motion MRI in healthy volunteers and patients with renal pathologies

To assess the feasibility of measuring tubular and vascular signal fractions in the human kidney using nonnegative least‐square (NNLS) analysis of intravoxel incoherent motion data collected in healthy volunteers and patients with renal pathologies.


| INTRODUCTION
DWI is a noninvasive MR technique for measuring the random motion of water on a molecular level and has previously been shown to provide differentiated information on renal microstructure and function. 1,2 Over the past few years, the mono-exponential ADC has emerged as a robust biomarker for various renal disorders such as renal artery stenosis, 3 renal fibrosis, 2 renal cell carcinoma (RCC), 4 ureteral obstruction, 5 renal failure, 6 and kidney graft rejection. 7 More recent studies, however, have reported that in intrinsically well-perfused organs (e.g., kidney, liver), the diffusion-attenuated MR signal is better described by the bicompartmental intravoxel incoherent motion (IVIM) model, which distinguishes between capillary perfusion and tissue diffusion. [8][9][10][11][12][13] Further, by incorporating a third component in the IVIM model, additional information about the kidney structure and filter function may be obtained. 14 Nevertheless, triexponential IVIM modelling in the kidney remains challenging due to usually limited b value sampling and low SNR of the diffusion-weighted images. 15 In the conventional biexponential IVIM, the pseudo-diffusion coefficient (D*) is often fixed [16][17][18] to improve the fit stability. However, fixing the intermediate (D* interm ) and/or fast (D* fast ) pseudo-diffusion coefficients when using the triexponential IVIM model may introduce bias to the perfusion fractions (f slow , f interm , f fast ). 15 As a promising alternative to the conventional approaches based on the nonlinear least-squares algorithm, a model-free nonnegative least-squares (NNLS) 19 method has recently been proposed for fitting the IVIM model in the brain, bone marrow, and kidney. 20,21,17,22 The NNLS analysis generates a spectrum with multiple peaks that represent diffusion and pseudodiffusion components with different diffusion coefficients and area fractions. 23 A key feature of this technique is that it does not require any a priori knowledge on input parameters and that the number of (pseudo-) diffusion components does not need to be specified in advance. 24 This might be particularly advantageous when studying biological tissues with a complex and heterogeneous structure, where the range of physiologically reasonable initial guesses is often unknown.
Previous studies confirmed the existence of 3 distinct diffusion components in the healthy kidneys. It has also been speculated that the intermediate-diffusing component originates from tubular structures associated with the pre-urine flow processes. 14,15 The volume fraction of tubules may change as a result of alterations in water filtration and reabsorption, tubular flow, or hyperglycemia. [25][26][27] From a clinical perspective, tubular atrophy and interstitial fibrosis, associated with decreased tubular volume fraction and increased interstitial volume fraction, are important histological features of progressive renal diseases. [28][29][30][31] Nevertheless, percutaneous kidney biopsy, which is currently considered the gold standard for diagnosis of nearly all pathologic kidney diseases, is invasive and can under-or overestimate the true severity of tubular atrophy and interstitial fibrosis for the entire parenchyma. 32 In contrast, IVIM MR imaging is a noninvasive method for assessing pathologic changes involving the whole kidney and could potentially become a valuable complement to the biopsy.
In this study, we applied the NNLS-based IVIM analysis in the human kidney (1) to assess the number of distinguishable diffusion components, and (2) to determine the renal tubular and vascular volume fractions in healthy volunteers and showcase patients with different kidney pathologies.

| Study population
The study was approved by the local ethics committee. After giving written informed consent, 12 healthy volunteers (5 females and 7 males; age range: 23-31 years; mean age 25.8 ± 2.7 years) with no history of renal diseases were included. The test-retest repeatability was assessed based on the data obtained in 5 subjects scanned twice with a time interval of 15 min and repositioning in the MRI. Furthermore, a patient with fibrotic kidney disease (male, age: 77 years, glomerular filtration rate: < 10 mL/min), a patient with a renal solid mass (male, age: 74 years, glomerular filtration rate: 71 mL/min), and a patient with a failed kidney transplant (male, age: 77 years, glomerular filtration rate: <10 mL/min) were examined. Histological examination of the solid renal mass was performed after radical nephrectomy. Subjects were not given any restrictions regarding fluid or food intake. 33

| In vivo MRI experiments
All MRI scans were performed on a 3 Tesla whole-body MRI system (Magnetom Prisma, Siemens Healthcare, Erlangen, Germany) using an 18-channel torso array coil and a 32-channel spine coil. To acquire anatomical images, a half-Fourier single-shot turbo spin echo sequence in all 3 image axes (axial, coronal, and sagittal) was used. Diffusionweighted images were obtained using a prototypical singleshot EPI sequence with reduced volume excitation (ZOOMit, Siemens Healthcare) 34

| Model and fitting
A measured diffusion signal y i can be described as a linear combination of multiple exponential decays 21 : where M is the number of diffusion coefficients D (in mm 2 /s), s j the relative amplitude for each D j , b i the b value (in s/mm 2 ), and N the total number of data points. The nonlinear optimization methods (e.g., Levenberg-Marquardt [LM] algorithm) solve Equation (1) for amplitudes s j at unknown D j and usually require an initial estimate of the solution to begin iterations. In contrast, the NNLS method is a noniterative approach and does not require any initial guesses. Whereas the nonlinear techniques often employ the restrictive assumptions of only a few exponential components, the linear inverse methods such as the NNLS approach assume a large number of known D j (usually between 100 and 200) and solve Equation (1) for the corresponding amplitudes s j , some of which might be 0. 23 In general, Equation (1) can be rewritten as follows 23 : where A ij is a constraint matrix representing the exponential decay functions. Because of noise contaminating y i , amplitudes s j cannot be exactly determined by simply taking the inverse of A ij . The misfit between the measured and modeled data can be quantified using the χ 2 statistic. In this study, the NNLS algorithm of Lawson and Hanson was employed to minimize the least-squares χ 2 misfit 23 : where all s j are implicitly constrained to be nonnegative s j ≥ 0. In output, NNLS yields a set of M delta functions with amplitudes s j for each logarithmically spaced D j . In order to produce a continuous distribution of diffusion coefficients S(D), which better represents a real physical system, and to improve fit robustness, a regularization term µ that minimizes distribution curvature can be incorporated into Equation (3) 35 : As a result, the NNLS fitting procedure provides a number of peaks representing the (pseudo-) diffusion components found in a given diffusion decay curve. The signal fraction is defined as the sum of S(D) in a specific range (D min to D max ) normalized by the sum of all S(D), whereas the mean diffusivity for this range is calculated as 35 :

| Data analysis
The trace-weighted DWI datasets were processed using in-house developed software (stroketool, Digital Image Solutions, Frechen, Germany). Retrospective motion correction was not performed in our study because it did not substantially improve the quality of our data. Only a single image slice containing the center of the right kidney was selected for further analysis because (1) it was not affected by respiratory motion compared with the remaining slices; (2) it contained the maximum parenchymal area, which allowed the largest cortical and medullary regions of interest (ROIs) to be obtained; and (3) both cortex and medulla were clearly visible thus enabling a reliable segmentation. Furthermore, the midslice of the kidney was assumed to be representative enough of the whole kidney in several previous DWI studies. [36][37][38] The datasets were smoothed with an isotropic 3 × 3 Gaussian kernel to increase SNR. Cortical and medullary ROIs and the tumor and cyst volumes were then manually segmented on the b 0 images by an experienced radiologist (a.l., 11 years of experience) in the open-source image segmentation toolbox ITK-SNAP (Version 3.8.0). 39 The measured diffusion-related signal decay was subsequently fitted using the NNLS method. Because multi-exponential fitting of noisy data are an ill-posed problem, a regularization parameter µ of 0.02 has been used as an input to make the NNLS fit more resilient. 35 Because NNLS algorithms perform better as SNR increases, 27 average signals from the cortex, medulla, and lesion ROIs were computed first and then analyzed using the NNLS technique. Mean and SD of the fitted diffusion parameters among the 12 healthy volunteers were obtained. In order to allow a direct visual comparison, voxel-wise NNLS fitting was also performed.
Furthermore, a triexponential IVIM model was fitted with the LM algorithm using the MatLab Curve Fitting Toolbox (MatLab R2018a, Mathworks). The initial parameter values used for the nonlinear fit are given in Supporting Information Table S2.

| Statistical analysis
The goodness of the NNLS fit was assessed with the coefficient of determination R-squared (R 2 ). Differences between the (1) cortical and medullary diffusion parameters were compared using Wilcoxon rank sum test. To examine the test-retest repeatability, within-individual coefficients of variations (CoV w ) were calculated as follows 40 : where x 1 and x 2 are the ROI-based values for the first and the second scan, ‼ x the mean value over all subjects, and N and n the total number of volunteers scanned twice and subject number, respectively.

| Subjects
All subjects were successfully scanned, and all datasets had sufficient data quality for further analysis. The solid lesion (maximum diameter: 7.3 cm) was classified by a pathologist as grade 4 clear cell renal cell carcinoma based on the 2016 World Health Organization/International Society of Urological Pathology classification. Histopathological images and a photograph of the resected kidney are displayed in Supporting Information Figure S1. Examples of the diffusion-weighted images at different b values and the ROI placement are shown in Figures 1 and 2, respectively.

| Parameter maps
Signal fractions of the diffusion components of an exemplary healthy kidney, the fibrotic kidney with multiple cysts, failed renal graft, and kidney with RCC and a large cyst are shown in Figure 3. In all investigated kidneys, the slow diffusion component had the highest percentage contribution. At the same time, the fraction of the slow component was much higher in the fibrotic cortex, failed renal graft, and RCC compared with the unaffected kidney. The f interm maps showed noticeable differences between healthy renal tissue and different kidney pathologies as well. The highest fractions of f interm were measured in the cysts that are filled with free-moving water. In contrast, the fibrotic cortex tissue, RCC, and renal graft parenchyma demonstrated decreased f interm contribution. The fast component appeared higher in high blood flow areas, for example, large blood vessels (see renal hilum of the fibrotic kidney). No contribution from f fast was found in the cysts and fibrotic parenchyma.

| ROI-based parameter analysis
The parameter values (and mean and SD, if possible) of the signal fractions and mean diffusivities are given in Table 1. The results of the triexponential IVIM fitting using the LM algorithm are presented in Supporting Information Table S3.
In general, the ROI-based results largely confirmed the findings of the voxel-wise analysis. The NNLS analysis revealed 3 distinct diffusion components in both renal cortex and medulla of the healthy kidneys. Signal fractions of the diffusion components declined with an increase in the corresponding mean diffusivity, and there were no statistically significant differences between the cortical and medullary diffusion parameters obtained for each component in the healthy kidneys (P > .05).
The cortical f interm value in the fibrotic kidney seemed to be slightly lower than that in the healthy renal cortex. The opposite was observed for the corresponding f slow values. Note that the fast decaying component f fast was not detected either in the fibrotic cortex or medulla. Compared with the healthy renal cortex, the MD interm value in the fibrotic tissue appeared higher. In the failed renal transplant, cortical and medullary f interm values were much lower than the corresponding values in the healthy kidneys. The f slow component was accordingly higher. At the same time, increased MD slow values were measured in the failed kidney graft compared with those obtained in the normally working kidneys.
The renal cysts were recognizable by substantially higher f interm and lower f slow than all other tissue types. As expected, no f fast component was measured. Both MD interm and MD slow values were lower than the mean diffusivities obtained in the healthy cortex and medulla.
The renal clear cell carcinoma had the highest f slow and lowest f interm among all investigated tissues. The MD fast values were also higher than those measured in the healthy kidneys.
Diffusion-weighted signal decays and corresponding NNLS fits obtained in the renal medulla of the healthy kidneys, fibrotic kidney, failed renal graft, cyst, and RCC are shown in Figure 4. The signal decay curve of healthy medulla represents data averaged over all volunteers. The coefficients of determination R 2 in all subjects were above 0.99, indicating high quality of the NNLS fits. The signal decay curves and NNLS spectra measured in the renal cortex are displayed in Supporting Information Figure S2.

| Test-retest repeatability
An example of signal fraction maps obtained by repeated scans on a healthy subject are shown in Figure 5. Visual inspection permits confirmation of the repeatability of the NNLS method for determining signal fraction of the diffusion components in the kidney. The differences resulting from a slight scan mismatch can be observed in the difference map obtained by subtracting 2 corresponding b 0 images. Mean values, SDs, and within-individual coefficients of variations (CoV w ) calculated in the healthy cortex and medulla for each diffusion component are given in Table 2. The highest repeatability of the signal fraction was achieved for the slow component, followed by the intermediate and fast components. Similarly, the repeatability of mean diffusivity increased with a decrease in diffusion coefficient.

| DISCUSSION
In the present study, the number and relative contributions of distinguishable diffusion components in healthy and pathological human renal tissues were assessed using the unbiased model-free NNLS analysis. In general, the NNLS analysis of a diffusion signal finds a set of discrete exponential components with amplitudes s j determined at diffusion coefficients D j . In the healthy kidneys, 3 distinct peaks at D j of about 10 -3 , 10 -2 , and 10 -1 mm 2 /s were detected in the NNLS spectrum and designated as fast, intermediate, and slow (pseudo-) diffusion component, respectively. Although the estimated values of the signal fraction and mean diffusivity were basically in line with those obtained in prior human and animal studies, 14,27 direct comparison is hindered primarily because of differences in b value distributions used at different sites. [41][42][43] According to the conventional IVIM model, the slow-and fast-decaying components are related to the molecular diffusion ("true diffusion") and tissue microcapillary perfusion ("pseudo-diffusion"). 44 The third component with diffusion rate being in the order of magnitude of free water has previously been reported in the liver, 45,46 prostate, 47 and kidney 14,15,21,27 and is generally believed to reflect free diffusion or microperfusion. Nevertheless, the correct interpretation of the observed diffusion components in the kidney IVIM still remains an open issue. A plausible explanation of their origin has recently been advanced by van der Bel et al. 15 They hypothesized that the intermediate compartment contains both the capillary perfusion and tubular pre-urine flow, which are in the same order of magnitude and therefore cannot be distinguished within 1 voxel. Moreover, the fast component, was assumed to originate from blood flow in larger vessels. Moreover, their data seemed to confirm previous findings that the slow compartment is most likely related to the passive water diffusion in the kidney parenchyma. 14,48 Our results tend to support this interpretation in terms of the measured diffusion parameters and their alterations under various pathological conditions. In the healthy kidneys, high values of f interm were observed in the renal cortex and renal columns (between the renal pyramids), which are the most vascularized parts of the kidney. 49 At the same time, the f slow map showed a pattern reflecting medullary pyramids with generally higher values than the surrounding tissue. The NNLS spectra obtained in the pathological renal tissues varied from those acquired in the healthy kidneys, most likely due to the differences in tissue cellularity and perfusion. Interestingly, the NNLS analysis resulted in less than 3 components being fitted in some cases. In fact, only 2 peaks at diffusion coefficients of about 10 -3 and 10 -2 mm 2 /s corresponding, respectively, to the slow and intermediate components were detected in the cyst and parenchyma of the fibrotic kidney. In the remaining tissues, high values of the third component at D j > 10 -1 mm 2 /s, indicating a very rapid decay of the diffusion signal, probably due to blood flow-related dephasing in larger renal vessels, were measured. 14,15 Indeed, an exemplary signal fraction map obtained for the fibrotic kidney showed a relatively high contribution of the fast component within the renal hilum. At the same time, both the f fast and f interm maps of the fibrotic parenchyma reflected capillary loss and hypoperfusion, which are the primary physical changes occurring in the fibrotic kidney. 50,51 A similar pattern was found in the failed transplant. This is consistent with previous studies, which showed a reduced perfusion fraction in the pathologic kidney transplant caused by lessened capillary vasculature and narrowed renal tubules. 52,53 In the cyst, there was no contribution from f fast but large f interm values were observed. Moreover, the estimated MD interm value of 3.39 × 10 −3 mm 2 /s was close to the self-diffusion coefficient of free water, which is around 3.0 × 10 −3 mm 2 /s at normal body temperature. 54 These findings are in line with earlier studies, where the intermediate pseudo-diffusion was speculated to be attributable to the freely diffusing water. 14,15,55 The investigated renal cell carcinoma demonstrated quite a different diffusion spectrum pattern compared to the healthy kidney. We argue that the decreased f interm and, at the same time, increased f slow , might be due to a mechanical hindrance of water motion caused by the diffusion-restricting elements such as cell membranes and macromolecules. 55,56 T A B L E 1 Diffusion parameters obtained from the NNLS-based IVIM analysis for each of the investigated tissue types.  In our preliminary study, unambiguous discrimination between different renal tissue types based only on a single diffusion parameter was not possible because signal fractions (and mean diffusivities) partially overlapped. This observation is in agreement with the literature. 56 We believe that a combination of diffusion parameters could potentially be used to reliably distinguish between different types of renal tissues/pathologies in the future. However, further studies are needed to explore the sensitivity of NNLS to changes in tubular and vascular flow and tissue cellularity. 27 There are several advantages of our study when compared with previous IVIM studies. First, our IVIM protocol allowed multi-b value DWI to be acquired within a clinically feasible time. To further reduce the scan time, an optimal acquisition strategy could be determined using simulations to select the most informative combinations of b values, gradient profiles, diffusion directions, and diffusion times. 57 Second, the NNLS method was able to identify a variable number of decay components without specifying their number beforehand. It has recently been speculated that when fewer components are present within a voxel than assumed in a model, overfitting can affect the parameter estimates. 22 This might be particularly relevant in pathological conditions, as shown by our preliminary results. However, an extensive comparison between the NNLS method and nonlinear fitting with mono-, bi-, and triexponential models, combined with histopathological analysis of different renal tissues, is needed to verify this hypothesis. Third, the NNLS approach does not generally require initial parameter guesses or a stopping criterion, which is probably its most important advantage over the most commonly used nonlinear least-squares based LM fitting. Good initial guesses and reasonable constraints  are crucial for avoiding convergence to local minima when fitting renal DWI data using the LM algorithm. Although the signal fractions and MDs derived from the triexponential IVIM model tended to vary slightly depending on the starting values used, they were generally in the same range of values as those obtained by the NNLS analysis (please see Supporting Information Table S3). Although the typical ranges of diffusion parameters in the healthy renal tissue are generally known, this might not always be the case when investigating different kidney pathologies. Lastly, we acquired IVIM data in patients with different kidney pathologies to directly compare spectral diffusion patterns obtained in various types of pathological and healthy renal tissues. Some study limitations need to be considered when interpreting our findings. First, the number of subjects was too low to permit a firm conclusion regarding the characterization of different renal pathologies. The potential of the NNLS method was shown only in 3 exemplary pathologies. In addition, the absence of statistically significant differences between the healthy renal cortex and medulla might be due to the low number of subjects and small effect size. Second, only a single central slice of the right kidney was analyzed in the current study. Assessment of the whole organ could generally provide more representative information on diffusion and perfusion. Moreover, the full coverage of a whole volume of interest might be particularly important for the evaluation of large and heterogenous renal tumors. Third, it was not possible to confirm the correlation between the diffusion parameters and the histopathological features of the fibrotic kidney and failed renal graft due to the lack of histological assessment. Fourth, although the repeatability of f slow , MD slow , f interm , and MD interm was satisfactory, relatively high coefficients of variation were obtained for the f fast and MD fast values. This is most likely due to the relatively sparse sampling and larger variations of the signal intensities in the lower b value regime. In fact, dense sampling of very low b values is required to precisely model the very fast initial signal decay. 44 Furthermore, despite using navigator-triggered prospective acquisition correction, 58 residual respiratory motion artifacts were still observed in some subjects. Hence, to reduce the motion-related variability of the observed components, only slices without any visible motion were analyzed. We decided not to perform any post-acquisition motion correction because it did not improve quality of our test datasets but could potentially introduce errors. As a result, the difference maps obtained by repeated acquisitions were slightly affected by a scan mismatch due to both the residual motion and respiratory phase mismatch. Last but not least, no dietary restrictions were imposed on the participants of our study, which might have contributed to the interindividual variations of the diffusion parameters. 59 All the above-mentioned issues should be addressed in future studies.

| CONCLUSION
We applied the model-free NNLS approach for spectral diffusion analysis of renal IVIM in the healthy and pathological kidneys. As expected, we identified 3 distinguishable components in the kidney DW signal, which may be related to 1) slow tissue diffusion, 2) intermediate tubular fluid flow and capillary blood flow, and 3) fast blood flow in larger vessels. Furthermore, our results indicate that the NNLS method is able to track changes in vascular and tubular volume fractions associated with various kidney pathologies. We believe that the NNLS analysis applied to renal IVIM data has the potential to become a valuable noninvasive diagnostic tool for measuring microscopic changes in renal tissues. However, further longitudinal studies and developments are required to explore the precise nature of the kidney diffusion signal components detected by the NNLS model.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section. Signal decay curves (left) and diffusion coefficient spectra (right) obtained in the cortex of healthy kidneys, a fibrotic kidney, and a failed kidney transplant. The NNLS spectrum shows distribution of detected exponential components, where each peak represents a (pseudo-) diffusion component with a mean diffusion coefficient (MD = geometric mean of peak) and signal fraction (f = area under the curve). The peaks at diffusion coefficient in the order about 10 -3 , 10 -2 , and 10 -1 mm 2 /s correspond to slow, intermediate and fast component, respectively TABLE S1 A summary of the IVIM acquisition details TABLE S2 Starting values and constraints used for the tri-exponential Levenberg-Marquardt fitting. The first parameter set was previously employed in the study of van der Bel et al, 15 and the second one uses the output parameters of the NNLS analysis as initial guesses TABLE S3 Signal fractions and mean diffusivity values obtained from tri-exponential fitting and non-negative least square (NNLS) analysis in twelve healthy volunteers and patients with various kidney pathologies. The tri-exponential fitting was performed using the Levenberg-Marquardt (LM) non-linear least squares algorithm with 2 different sets of starting values and constraints as defined in Supporting Information