Theory

In view of our aim to measure medullary perfusion, i.e., blood flow from cortex to medulla, we identify the simplest model that can correctly describe the exchange of contrast agent between cortex and medulla (Fig. 1). A contrast agent molecule can follow essentially two trajectories between these two regions:

(1)

molecules that are filtered out of the plasma compartment (PA, cortex) enter the proximal tubules (PT, cortex) then the loop of Henle (LH, medulla), the distal tubules (DT, cortex) and the collecting ducts (CD, medulla) before being evacuated as urine (U);

(2)

molecules that are not filtered out enter the vasa recta (VR, medulla) then the venous plasma (PV, cortex) before being evacuated as venous blood (V), or else they stay in the cortex and go straight to the venous plasma (PV, cortex).

Fig. 1

Schematic representation of the seven-compartment model. Flows with tracer are denoted FX,Y where the flow is from compartment X to compartment Y or by an extraction fraction times Fcor (the arterial flow into the cortex). Resorption or tracer free flows are expressed as fractions fX of FT (EFFFcor), with X the compartment from which resorption occurs. Note that the arterial plasma compartment is not strictly arterial since it also includes glomerular and peritubular capillaries. A renal arteries, PA arterial plasma compartment, PV venous plasma compartment, V renal veins, VR vasa recta compartment, PT proximal tubules compartment, LH loop of Henle compartment, DT distal tubules compartment, CD collecting duct compartment, U urine/renal calyces

We end up with a 7-compartment model (Fig. 1) with 3 plasma compartments (PA, VR, PV) and 4 tubular compartments (PT, LH, DT, CD). The model is fully defined by 10 free parameters, summarized in Table 1. They are the mean transit times (MTT) of each of the 7 compartments, the cortical perfusion Fcor, the medullary extraction fraction Emed, and the glomerular extraction fraction EFF. The target parameter medullary perfusion Fmed can be derived as Fmed = Emed(1-EFF)Fcor.

Table 1 A list of primary and derived model parameters with units and their description. MTT: mean transit time.

If we write MTTX for the MTT of a compartment X, then we can define the propagator HX and residue functions RX of the compartment as follows:

$$R_ = e^}} H_ = \frac}}}}$$

With these definitions, the solutions of the model equations can be written out compactly by following the trajectory of the tracer flux from the arterial inlet to the compartment. For a given arterial concentration CA(t), the concentration in each of the 7 compartments is given by (\(*\) denotes convolution):

$$C_ \left( t \right) = \frac }}R_ *F_ C_ \left( t \right)$$

(1)

$$C_ \left( t \right) = \frac }}\left[ *H_ *\left( } \right)\left( } \right)F_ C_ \left( t \right) + R_ *H_ *H_ *E_ \left( } \right)F_ C_ \left( t \right)} \right]$$

(2)

$$C_ \left( t \right) = \frac }}R_ *H_ *E_ \left( } \right)F_ C_ \left( t \right)$$

(3)

$$C_ \left( t \right) = \frac }}R_ *H_ *E_ F_ C_ \left( t \right)$$

(4)

$$C_ \left( t \right) = \frac }}R_ *H_ *H_ *E_ F_ C_ \left( t \right)$$

(5)

$$C_ \left( t \right) = \frac }}R_ *H_ *H_ *H_ *H_ *E_ F_ C_ \left( t \right)$$

(6)

$$C_ \left( t \right) = \frac }}R_ *H_ *H_ *H_ *E_ F_ C_ \left( t \right)$$

(7)

It is unlikely that a model of this complexity can be identified based only on a measured curve over the entire kidney parenchyma. Instead, we follow a strategy introduced by Baumann et al. [10] of collecting concentrations in the cortex and medulla separately. Since medullary perfusion quantifies the exchange between these two regions, we hypothesized that this would produce a well-defined inverse problem. The cortex consists of the arterial plasma, venous plasma, proximal tubules and distal tubules:

$$C_ \left( t \right) = V_C_ \left( t \right) + V_C_ \left( t \right) + V_C_ \left( t \right) + V_C_ \left( t \right)$$

(8)

The (fractional) volumina correct for the relative contribution of each compartment and cancel out in the final equation. The medulla consists of the vasa recta, loop of Henle, and collecting ducts:

$$C_ \left( t \right) = V_C_ \left( t \right) + V_C_ \left( t \right) + V_C_ \left( t \right)$$

(9)

Patient data collection

Patient data were selected from an ongoing study (iBEAt [11]) on imaging biomarkers in diabetic kidney disease. The main inclusion criteria were type-2 diabetes, eGFR > 30 ml/min/1.73m2. All subjects signed informed consent prior to inclusion. For the purposes of this paper, the first 24 patients of one recruiting site were selected, excluding cases where DCE-MRI data were deemed to be of insufficient quality. Quality control consisted of visual inspecting of the raw images in a DICOM viewer. Any images that show artifacts, signal drop-outs, data truncation or extreme noise levels suggesting coil issues, were excluded.

Patients were scanned on a 3 T MR system (Magnetom Prisma, Siemens, Erlangen, Germany), feet first with an 18-channel phased array body coil. Apart from DCE-MRI, the scan protocol included a range of sequences including DIXON, T1/T2/T2*-mapping, diffusion weighted imaging and diffusion tensor imaging, ASL and phase contrast. DCE-MRI was performed in free breathing with a 2D fast gradient echo sequence with a non-selective saturation pulse prior to each slice readout. 2D acquisition was chosen because of the smaller acquisition time per slice as compared to a 3D volume acquisition. This renders the images less vulnerable to (respiratory) motion artifacts. Eight oblique slices through the kidney and one transverse slice through the descending aorta were acquired with acquired voxel size of 2.78 × 2.08 × 7.5 mm3 reconstructed to a voxel size of 1 × 1 × 7.5 mm3. The field of view was 400 × 400mm and a parallel imaging factor of 2 was employed, along with a partial Fourier factor of 7/8. Temporal resolution was 1.6s and total acquisition time was 7 min. Flip angle was 10 degrees, TR and TE were 179 ms and 0.97 ms, respectively, with an echo spacing of 2.2 ms. Image acquisition started 20s before bolus injection of 0.05 ml/kg gadoterate meglumine at a rate of 2 mL/s followed by a 20 mL saline flush at a rate of 2 mL/s.

DCE-MRI images were aligned slice-by-slice using model-driven registration [12]. Registration accuracy was evaluated by visual comparison of the coregistered dynamics against the original time series, and by comparing motion corrected versus uncorrected datasets. Motion correction results were rejected if the coregistration had created unrealistic deformations or had blurred out detail that was visible in the source images. Segmentation was performed on the registered data using k-means clustering [13]. Three clusters were created: one to represent the cortical voxels, one for inner medullary voxels and one for voxels containing partial voluming/outer medulla (on average 22% of all voxels labeled as kidney parenchyma). The latter were not used in further processing. Segmentation of the aorta for calculation of the arterial input function was performed by k-means clustering as well. Any voxels outside of the selected anatomical regions were manually removed. Average cortical, inner medullary, and aorta time–intensity curves were calculated from each data series.

Implementation of the model

The model summarized in Eqs. 8 and 9 was implemented in Matlab (R2019a, MathWorks, Natick, MA, USA), using signals in arterial blood, cortex and inner medulla as input data and producing measurements of all 10 free parameters as outputs. Arterial and tissue concentrations were assumed to be proportional to the signal change with respect to the baseline level. A hematocrit of 0.45 was assumed to convert the aorta blood concentration to plasma concentration. For a given arterial plasma concentration, predicted tissue concentrations were calculated using Eqs. 8 and 9 with convolutions calculated as described in [14]. The model parameters were fitted to inner medullary and cortical signal curves simultaneously in a single regularized fit [15]. The model fit was repeated 100 times with different initial parameter values randomly chosen within physiological ranges, but without bounds on the fits. Of the 100 obtained solutions, only those solutions were selected with all parameters within physiological bounds. The solution with the lowest chi-square was selected as the result of the fit. The procedure is graphically depicted in Fig. 2.

Fig. 2

The fitting procedure depicted schematically. As part of the fitting procedure, a regularized unconstrained fit is performed 100 times with randomly chosen initial values. Next, the solutions within physiological bounds are selected. Last, the solution with the lowest chi-squared is selected

Simulated data collection

The aorta concentrations measured in the patients were used to generate kidney tissue signal–time curves with known ground-truth model parameters (Fig. 3). To determine realistic ground truth and bounds for all parameters, literature values were used where available. These were subsequently adjusted by data obtained by explorative model fits on patient data. Sets of ground truth values for all parameters were obtained from the normal distribution with a standard deviation of 20% around those theoretical values. Using those sets of ground truth values and the measured arterial plasma concentrations, signal–time curves were generated. Gaussian noise with a fixed standard deviation was added to the signals to obtain an arterial contrast-to-noise ratio (CNR) of 150, similar to the patient data. Arterial CNR was defined as the ratio between peak arterial signal change and baseline standard deviation.

Fig. 3

Data generation to determine stability and accuracy of the fitting procedure. Left: simulated data: (1) for each available arterial input function (aorta concentrations measured in patients) 4 sets of kidney tissue signal-time curves were generated, for a total of 96 unique curves with unique random choices of parameters, (2) the fitting procedure as depicted in Fig. 2 was applied 50 times to determine fitting stability and accuracy. Right: patient data: (1) for each patient, two datasets were available (left and right kidney), (2) the fitting procedure as depicted in Fig. 2 was applied 50 times to determine fitting stability

Data analysis

While the analysis is designed primarily to measure inner medullary perfusion, the 9 other parameters generated by the analysis are reported as well as they can be of use for additional quality control.

Patient data were used to verify that the model fitted the data well, compare results against known literature values, test sensitivity of the results to the choice of initial values, and compare shared parameters against a common reference model. Goodness of fit was assessed visually by inspection of fits overlaid on data. Results for all model parameters were reported as groupwise mean (standard deviation) and compared to literature values where available. In order to test model sensitivity to initial values, the fitting procedure was repeated 50 times for each dataset. This procedure is schematically depicted in Fig. 3. The coefficient of variation (CoV, defined as standard deviation/mean) was calculated over these 50 fits. If the iterated fit is perfectly stable, this CoV equals zero. Values for tubular flow FT (GFR per unit of tissue volume) and Fcor were compared to values obtained using the current reference two-compartment filtration model [7]. The intraclass correlation coefficient ICC was used to measure agreement between both measurements (two-way mixed effects, consistency, single measurements).

The simulated data were used to determine sensitivity to initial values and accuracy of the model fit if used with noisy data. The bias with respect to the ground-truth value was used to determine accuracy of the model fit. To determine whether the model is sensitive enough to detect individual differences in Fmed, fitted Fmed was compared to true Fmed for all ground truth curves using the intraclass correlation coefficient (ICC, two-way mixed effects, absolute agreement, single measurement). Sensitivity of all parameters to initial values was measured in a similar way as in patient data by repeating the fitting procedure 50 times for each simulated dataset (Fig. 3) and quantify sensitivity using a CoV defined as standard deviation/mean.