Quantitative transport mapping (QTM) of the kidney with an approximate microvascular network.

PURPOSE: Proof-of-concept study of mapping renal blood flow vector field according to the inverse solution to a mass transport model of time resolved tracer-labeled MRI data. THEORY AND METHODS: To determine tissue perfusion according to the underlying physics of spatiotemporal tracer concentration variation, the mass transport equation is integrated over a voxel with an approximate microvascular network for fitting time-resolved tracer imaging data. The inverse solution to the voxelized transport equation provides the blood flow vector field, which is referred to as quantitative transport mapping (QTM). A numerical microvascular network modeling the kidney with computational fluid dynamics reference was used to verify the accuracy of QTM and the current Kety's method that uses a global arterial input function. Multiple post-label delay arterial spin labeling (ASL) of the kidney on seven subjects was used to assess QTM in vivo feasibility.
RESULTS: Against the ground truth in the numerical model, the error in flow estimated by QTM (18.6%) was smaller than that in Kety's method (45.7%, 2.5-fold reduction). The in vivo kidney perfusion quantification by QTM (cortex: 443 ± 58 mL/100 g/min and medulla: 190 ± 90 mL/100 g/min) was in the range of that by Kety's method (482 ± 51 mL/100 g/min in the cortex and 242 ± 73 mL/100 g/min in the medulla), and QTM provided better flow homogeneity in the cortex region.
CONCLUSIONS: QTM flow velocity mapping is feasible from multi-delay ASL MRI data based on inverting the transport equation. In a numerical simulation, QTM with deconvolution in space and time provided more accurate perfusion quantification than Kety's method with deconvolution in time only.

INTRODUCTION

The traditional approach to quantitative perfusion mapping uses Kety’s equation, where the arterial input function (AIF) at each voxel is not measurable and is assumed to be a global function ignoring dilution, delay, and dispersion. For example, in arterial spin labeling (ASL) that uses radiofrequency (RF) labeled water as tracer, the AIF is assumed at the upstream arterial site where RF label is applied. The use of global AIFs intrinsically causes errors in quantitative perfusion mapping. Although there have been attempts at modeling and correcting for delays and dispersions, perfusion quantification using Kety’s equation remains prone to errors of local arterial dispersion and has never been validated against a ground truth.

Regional blood flow in tissue is fundamentally governed by the transport equation of mass and momentum fluxes, or the continuity equation and Navier‐Stokes equation. , The solution to the transport equations requires knowledge of the detailed vascular anatomy at a sub‐voxel scale. Integration over the voxel of the detailed transport solution can be used to model signals acquired in macroscopic time‐resolved imaging of tracer passage in tissue. The inverse solution to this transport model allows determination of transport parameters including the velocity vector field. , Accordingly, we are developing this quantitative transport mapping (QTM) to quantify tissue perfusion from time‐resolved imaging of tracers.

We report here the feasibility of QTM for the human kidney where the fractal microvasculature model has been established for the intravoxel blood vessel structure. Renal blood velocity/flow maps obtained using QTM in numerical simulation agree with the ground truth, and its in vivo feasibility was assessed in healthy volunteers using renal ASL MRI data with multiple post label delays (PLDs).

THEORY

Blood tracers, including the non‐invasively labeled blood water for ASL and injected contrast agent for dynamic contrast enhancement (DCE) and dynamic susceptibility contrast (DSC), are assumed to have the same transport parameters including the velocity vector field as blood. , Here, we review the transport equation for the labeled blood water in ASL as a forward problem and describe the corresponding inverse problem.

Forward problem

The forward model of QTM includes the microvascular network, Navier‐Stokes equation, and the continuity equation. Given the microvascular model, the intravascular blood velocity is governed by the steady‐state Navier‐Stokes equation and the tracer concentration by the continuity equation. where denotes the density of the blood (1.06 g/ml), the pressure, the blood viscosity, the domain of intravascular space, the tracer concentration, the time derivative, the spatial gradient, the blood velocity field, the position vector, the diffusion coefficients, and the tracer decay rate. with 1650 ms for labeled blood tracer data in ASL on 3T scanner. , Pulsation is ignored for simplicity. The velocity field is regarded as constant in time and the diffusion coefficient is assumed to have a negligible effect, that is, . The velocity in the tracer transport equation Equation (3) is assumed to be the same as that of blood flow in Equation 1.

The boundary conditions for Equations (1 and 2) are determined by the microvascular model and the blood pressure difference between the inlet (main arterial) and the outlet (capillary level vessel at arterial side). The boundary conditions for Equation (3) are determined by the tracer supply time duration and the permeability of the vessel wall. For simplicity, we assume the vessel wall is impermeable, that is, the tracer is confined within the intravascular space, although permeability and other processed can be included in the above equation. Equations ((1), (2), (3)) together serve as the forward model of the proposed QTM to provide the solution of blood velocity and tracer concentration.

In summary, tracer passage through tissue as captured in dynamic perfusion imaging is fundamentally governed by the physics of transport equation consisting of mass and momentum flux and boundary conditions of tissue microvasculature.

Inverse problem

The inverse problem in QTM is to reconstruct the voxel‐based blood velocity/flow using the tracer concentration image data at the voxel level of spatial resolution. Approximating in studying convection dominated transport phenomena, we solve for that satisfies Equation (3) for a given multiple time frames of tracer concentration . As the ASL tracer signal in MRI is weak, approximately 2 ~ 3% of the total signal in the brain and 5 ~ 7% in the kidney, and is prone to noise and motion induced errors, regularization is used to reduce their effects.

Let’s first consider the inverse problem in the continuous space. Given the measured 4D time‐resolved tracer data , the inverse problem of fitting to Equation (3) for estimating the velocity can be formulated as the following minimization problem.where the tracer concentration at the ‐th time frame , the number of time frames of concentration data used, the standard Euclidean norm, the norm, and the regularization parameter. can be determined using the L‐curve method to best reduce the ill‐posedness of the inverse problem. The blood flow within a blood vessel can be calculated from the reconstructed velocity field by integrating the velocity on the corresponding cross‐sectional area of the blood vesselswhere is the cross‐sectional area of the investigated blood vessel.

Now we introduce discretization at the imaging voxel level, assuming that the microvasculature is known. The notation for the voxel discrete space variables is used to distinguish them from the microvascular continuous space variables . Equation (3) needs to be integrated over a voxel space, where are the voxel indices, and are the voxel dimensions:

Here, denotes the volume of the voxel. Define a binary microvasculature mask if is in the intravascular space and otherwise. Then Equation (6) becomes.

Note that the microvasculature is assumed to be time‐independent, that is, .

Let us define the intrinsic concentration at voxel as and the intravascular blood volume. Applying the divergence theorem, Equation (7) becomes.where denotes the surface of the voxel space and is the outward unit normal vector of . With the assumption of smooth velocity field and concentration , Equation (8) can be rewritten as the following voxelized continuity equation (see Supporting Information documents, which are available online, for details):with defined as

The solution of Equation (9) is a quantity across the voxel surfaces with unit , according to the definition in Equation (10). Therefore, in the discrete space of imaging voxels, the transport equation after integration over a voxel (Equation 9) appears in an apparent form similar to that in the continuous space (Equation 3).

According to Supporting Information Equation (S9) (see Supporting Information documents for details), blood flow on each surface of the voxel is.

If we define the following vasculature parameter.where the superscripts denote the Right, Anterior, and Superior surfaces of a voxel, respectively (see Supporting Information Figure S1 for details). Then the average blood velocity on the voxel surfaces can be derived as

Analogous to Equation (4), QTM solves the blood flow velocity according to the inverse problem of Equation (9), representing deconvolution in both space and time:

The blood flow can be determined according to the velocity and voxel surface vasculature:

In summary, the voxelized transport equation can be used to fit voxel‐averaged concentration measured on space‐time resolved imaging data, and the fitting or the inversion of the voxelized transport equation can determine an apparent velocity vector field defined by averaging on the voxel surface.

Renal microvascular network construction

A renal microvasculature for each subject was constructed in the following way. The 3D shape of the kidney was obtained from a segmentation of the mean ASL data in time. The MRA was used to segment the renal arteries using ITK‐SNAP semi‐automatically and to identify the directions of their blood flow. Starting from the segmented renal artery branches, a microvasculature structure in the segmented 3D kidney shape was generated using a bifurcation branch rule. The iterative process to generate renal microvascular branches was continued until the renal blood volume in 85% of the voxel within the cortex reached fell within the 16%‐20% range. , The output microvasculature will be used as the fluid space for the forward problem computation based on the fluid dynamics. The structural parameters, that is, the intersection area , and between vessel and voxel surface and the blood vessel volume , were calculated for all the voxels. These parameters will be used for both numerical simulation experiments and in vivo experiments.

QTM numerical simulation with microvascular network

The numerical simulation was done based on the simulated microvascular network obtained for one of the subjects recruited in this study (see above). In order to verify the proposed QTM method, a ground truth blood velocity and pressure in this microvascular network was computed following the method proposed by Lorthois et al. In this approach, the 3D microvascular network was simplified as a 1D non‐linear network of cylindrical model (see Supporting Information Figure S2 for details). The flow rate through each vessel segment was related to the pressure drop by Poiseuille’s law, which is a special case of the Navier‐Stokes equations Equations (1 and 2) (see Supporting Information document for details). Therefore, our simulations are based on the Poiseuille’s flow:withwhere denotes the pressure at node and is the hydrodynamic conductance of segment , is the apparent viscosity of the blood, and are the segment radius and length, respectively. Note that depends on the vessel diameter and hematocrit as shown in Lorthois et al. Applying Equation (16) to every segment of the microvascular network and considering the conservation of mass at each interior node , we have a large sparse linear system

Given the flow rate at each boundary node, Equation (18) can be solved using the conjugate gradient (CG) method to get the pressure at each interior node. The blood flow rate is calculated using Equation (16).

Consequently, the average velocity in each vessel segment can be calculated by.

The parabolic velocity profile of the segment is, therefore, expressed as.where is the radial distance from the center of the cylindrical vessel.

For the numerical simulation in the microvasculature, the blood pressure in Equation (18) was solved by giving a boundary condition of blood flow rate at the boundary of the microvascular network, that is, the flow inlet and outlets. The inlet flow was estimated by multiplying artery blood velocity (taken from literature measurements ) with the artery cross‐sectional area estimated from MRA. The outlet flow was assumed to equal for all the terminal capillary segments and the sum of all outlet flows was set equal to the inlet flow. Then the blood flow and velocity for each vessel segment were computed using the calculated pressure. The ground truth of the voxelized flow rate was calculated by averaging the inflow and outflow for each voxel.

The tracer concentration distribution in the microvasculature can be computed by solving the partial differential equation (PDE) Equation (3) with the velocity in Equation (20) for multiple time steps. However, the approach of solving Equation (3) using the finite element method (FEM) or finite difference method (FDM) in the large vascular system is challenging and time‐consuming. We adopted the approach of the analytical solutions by Chen et al and Massabo et al , to compute the concentration distribution in the microvascular network (see Supporting Information Figures S2 in the Supporting Information and Figure A1 in the Appendix for details). The approach of the analytical solution is summarized in the Appendix.

FIGURE 1

Illustration of microvasculature construction and structural parameters. A, Renal MRA image. B, Constructed microvasculature in the kidney using the renal arterial vessel information and shape information from MRA image in (A). C, Calculated renal blood vessel volume map (unit:) from constructed microvasculature in (B), D‐F, Calculated intersectional areas(unit:) of renal blood vessel on the right, anterior, and superior surfaces of a voxel, respectively

METHODS

The proposed QTM was applied to 4D time‐resolved ASL MRI data. All subjects satisfied the MR safety criteria. Informed consent was obtained from each subject under approval from the institutional review board at our institution.

Data acquisition

Pseudo‐continuous ASL (PCASL) 3D fast spin echo (FSE) data were acquired in the kidney in healthy subjects (N = 7; age 27 ± 6, six male, one female) using a GE MR750 3T scanner (GEHC, Milwaukee, WI) using a 32‐channel body coil. Acquisition parameters include 2.5 × 2.5 × 4 mm3 voxel size, 128 × 128 × 36 matrix, 10.5 ms echo time, 111° flip angle, three signal averages, ~4.5 min scan time. Four PLDs (PLD = 1025 ms, 1525 ms, 2025 ms, 2525 ms) were used during the data acquisition to obtain temporal resolution. Background suppression and synchronized breathing were used to enhance the signal quality and reduce motion artifacts. , Next, a non‐contrast 3D MR angiogram (MRA) data were acquired for each subject, using acquisition parameters: IFIR sequences, 0.625 × 0.625 × 2 mm3 voxel size, 512 × 512 × 128 matrix, 4 ms repetition time, 2 ms echo time, 50° flip angle, and respiratory gating.

Image analysis

QTM was performed by solving Equation (14) for blood flow vector field and Equation (15) for absolute blood flow on all data. Rigid image registration was performed for all PLDs using the first frame as reference. The intrinsic concentration in Equation (14) was derived from the motion‐corrected 4D ASL data as follows.where is the renal blood volume. Here, we assumed a linear relationship between image intensity and tracer concentration. The optimization problem Equation (14) was solved using the alternating direction method of multipliers (ADMM) with the conjugate gradient algorithm as a subroutine for solving the linear sub‐problems. The regularization parameter was chosen by performing an L‐curve analysis.

The solution of Equation (14) is the flow passing through three surfaces of a voxel. The blood flow map from QTM can be obtained directly from Equation (15)in (mL/100 g/min), where denotes the scaled flow of Equation (15), and is the density of the tissue.

The blood flow () map (mL/100 g/min) was computed by fitting ASL multiple delay data to the following signal model , , where is the difference between control image and labeling image at , the proton density‐weighted image, the labeling efficiency, the tissue/blood partition coefficient, , in renal cortex and in the medulla of kidney, the arterial transit time (ATT) of blood from labeling location to the kidneys estimated from the weighted delay defined in Wang et al and Dai et al, , the labeling time. Using Equation (23), the AIF for all the voxels was assumed to be the arterial concentration at the aorta site of spin labeling corrected by the decay. As suggested Koh et al, to reduce noise and enhance the signal‐to‐noise ratio (SNR) of the image, the data were registered to the to reduce motion artifacts and smoothed using a smoothing filter with a diameter of 8 mm. Masks for the renal cortex and medulla regions were segmented for each subject using thresholding the average ASL data across PLDs. The blood flow was obtained from deconvolution in time by performing nonlinear least square fitting of Equations (23) for all PLDs. Average blood flow in the renal cortex and medulla was computed for both QTM and Kety’s method. For reference, Equation (23) for Kety’s flow as derived by Buxton is a time convolution of the differential Kety’s equation modified with a decay term:where AIF(t) is , , and as in Equation (3).

RESULTS

Figure 1 shows the microvasculature constructed from MRA and M0 data and the corresponding structural parameters calculated from the microvascular network using voxelization. Figure 1A shows the MRA image of the kidney, Figure 1B the constructed microvascular network, Figure 1C the absolute blood volume, and Figure 1D‐F the vessel intersectional areas on the voxel surfaces.

Figure 2 demonstrates the simulated pressure, velocity, flow in the numerically generated microvascular network: Figure 2A shows the pressure in mmHg computed using Equation (18), Figure 2B the velocity in mm/s computed using Equation (20), Figure 2C the flow in mL/min computed using Equation (16). These results were used to generate the ground truth of velocity and flow maps at the voxel level. Figure 2D‐G presents the simulated tracer concentration at various time points with tracer labeling duration 1.5 s. The parabolic velocity field in Figure 2B was used for the tracer delivery simulation.

FIGURE 2

Simulated pressure in Equation18(A), velocity in Equation20(B), flow in Equation16(C) in the constructed microvascular network of a kidney. D‐G, Simulated concentrations at four time frames in the simulated microvascular network of a kidney using the velocity field from Equation20

Figure 3 shows the voxelized simulated concentration data at 4 PLD as shown in Figure 2D‐G.

FIGURE 3

A‐D, Voxelized simulated concentration data of Figure 2D‐Gat the voxel size of [4, 2, 2], which is typical ASL data voxel size. This voxelized data will be used for the QTM simuation reconstruction

Figure 4 shows the comparison of blood velocity and blood flow between ground truth, QTM, and Kety’s method in the numerical simulation. Figure 4’s top panels (A), (B), and (C) show the ground truth velocity magnitude, QTM reconstructed velocity magnitude, and QTM velocity magnitude error, respectively. The bottom panels (D‐H) show the ground truth blood flow, QTM reconstructed blood flow, QTM blood flow error, Kety’s method blood flow, and Kety’s method blood flow error, respectively. The root mean square errors (RMSE) in the estimated flows were 49.7 mL/100 g/min and 121.8 mL/100 g/min, for QSM and Kety’s method, respectively. This corresponds to 18.6% and 45.7% of the mean of the ground truth flow 266.6 mL/100 g/min, respectively.

FIGURE 4

The numerical results comparison between QTM, Kety’s method and the ground truth. The top panel (A‐C) shows the ground truth velocity magnitude, QTM reconstructed velocity magnitude and QTM velocity magnitude error, respectively. The bottom panel (DH) show the ground truth blood flow, QTM reconstructed blood flow, QTM blood flow error, Kety’s method blood flow, and Kety’s method blood flow error, respectively

Figure 5 presents the in vivo ASL data for healthy subject at four PLDs.

FIGURE 5

In vivo ASL data on healthy subject 2. A‐D, The ASL data at four PLDs = 1025 ms, 1525 ms, 2025 ms, 2525 ms (left to right, top to bottom)

Figure 6 displays results in a healthy subject. Figure 6A shows a 3D view of the QTM velocity field in one slice. Figure 6B shows a comparison of the renal blood flow estimated using QTM (top) and Kety’s method (bottom). QTM provided a more uniform RBF in the renal cortex compared to Kety’s method. Across the seven healthy subjects, QTM provided mL/100 g/min and Kety’s method gives mL/100 g/min for the cortex region; QTM provided mL/100 g/min and Kety’s method mL/100 g/min for the medulla region.

FIGURE 6

QTM results on healthy subject. A, The 3D display of the vector field from QTM at one slice. B, The renal blood flow comparison between QTM (top) and Kety’s method (bottom)

Figure 7 shows the statistical analysis between the results from QTM and Kety’s methods. The top (bottom) row show the results for the cortex (medulla) region. Linear regression shows that coefficient in cortex region is 0.99 and in medulla region is 0.37. Bland‐Altman analysis shows that the results from QTM are underestimated compare to that of Kety’s method.

FIGURE 7

Statistical analysis between the results from QTM and Kety’s methods. The top (bottom) row show the results for the cortex (medulla) region. A, The correlation of the results between QTM and Kety’s method in the cortex region. B, The Bland‐Altman plot of cortex blood flow from QTM and Kety’s method. C,D, The corresponding plot of (A) and (B) in the medulla region

Table 1 shows the mean and SD of the measured blood flow in cortex and medulla using QTM and Kety’s methods for all the seven subjects. The coefficient of variation (CoV) in cortex region of all the subjects are also listed in the Table 1, CoV in cortex is 0.25 from QTM vs. 0.33 from Kety’s method.

TABLE 1

The regional blood flow results from QTM and Kety’s method, including blood flow in the cortex and medulla for seven healthy subjects and their mean and SDs

Method	Region	S1	S2	S3	S4	S5	S6	S7	Mean	SD
QTM	Cortex	443.78	518.33	445.46	391.19	417.43	368.55	518.84	443.37	58.19
	Cortex CoV	0.26	0.28	0.17	0.32	0.23	0.28	0.23	0.25
	Medulla	153.89	246.42	310.96	68.88	146.15	119.58	287.51	190.49	91.45
Kety	Cortex	486.16	553.02	471.73	442.70	459.52	415.83	546.37	482.19	51.22
	Cortex CoV	0.35	0.31	0.34	0.36	0.30	0.32	0.30	0.33
	Medulla	257.80	326.74	211.35	105.48	234.14	242.50	313.97	241.71	73.31

Cortex CoV is the coefficient of variation in the cortex region to measure the homogeneity of the blood flow distribution in the cortex.

DISCUSSION

Our numerical simulation and preliminary in vivo results demonstrate the feasibility of QTM to reconstruct the blood flow map and average velocity map from multiple delay ASL data in the kidney. The QTM approach is based on the biophysical description of blood transport in the tissue microvascular network. Voxelization with the microvasculature discretizes the continuous model for fitting time‐resolved imaging data. QTM is the solution to the inverse transport equation from imaging data to blood flow at the voxel level. Numerical results show that the reconstructed blood flow and velocity maps from QTM agree with the ground truth, and that the large error in the traditional Kety’s method is reduced by more than two‐fold using QTM.

Accurate mapping of regional blood flow in medical imaging is known to be very challenging. Although there are researches comparing ASL to positron emission tomography (PET) or microspheres for blood flow mapping, , it has not been validated against the ground truth blood flow directly since the ground truth in the real world is unknown. , , And the gold standard microsphere deposition method or PET method suffer from many assumptions including the use of a global AIF. The simulation in the current work based on the transport equation and microvascular network provides a ground truth for assessing quantification accuracy. The QTM approach does not assume an AIF as in the traditional Kety’s approach, , but requires detailed knowledge of intravoxel microvasculature, suffering consequently from the imperfect knowledge. The error in QTM flow may be caused by the voxelization error of concentration data, because the voxelization could possibly increase the error of spatial and temporal derivative of the tracer. This error could be reduced by increasing the data resolution and accordingly decreasing the voxel size. Errors in Kety’s method may come from the assumption of a single global AIF (see Figure A2 in the Appendix for details). In vivo results show that even though the results from both methods agree with literature values, the blood flow map obtained using QTM with the microvascular network approximation generates a more homogeneous flow expected in the renal cortex region compared to Kety’s method (mean CoV in cortex for QTM 0.25 vs. Kety 0.35 as shown in the Table 1).

QTM is the first method that uses the transport equation to quantify the blood flow at the image voxel level to the best of our knowledge. Although several groups have proposed to use the idea of mass conservation in PDE (Equation 3) to analyze the blood transport phenomenon, most of them focus on the forward problem analysis. , , , The 2D inverse problem of Equation (3) was considered in tumor assuming only (constant) diffusion. The assumption of the constant diffusion coefficient obviously deviates from the reality, and the direct inversion method for solving for the diffusion coefficient is likely prone to the noise in the data, especially with second derivatives of noisy data in the denominator. The direct inversion approach in Ref. 40 is not applicable to the blood perfusion case since blood velocity in 3D is directional. Our group has proposed to use porous medium approximating tumor tissue and solved the inverse problem to quantify the blood velocity using PDE Equation (3) in 3D. This work is an extension of the previous work by modeling an approximate microvascular network in the tissue to better quantify blood flow.

There are several benefits of the proposed QTM method. The QTM model is a biophysically based method derived from mass transport equation accompanied by rigorous mathematical derivation. Given the microvascular information, no AIF is required, as is the case with Kety’s method of quantitative perfusion mapping, which generalizes the Fick’s principle for an organ to mass conservation in a voxel with several assumptions. Kety’s method relies on empirically determined parameters, such as labeled blood arrival time, labeling efficiency, blood/tissue partition coefficient, labeling plane position, etc. , , A small change of these parameters can result in a large change in the blood flow value. Kety based methods have been proposed that do rely on many empirical parameters or assumed global AIF, but the fact remains that a single global AIF, even those reconstructed from blind convolution, are not appropriate for the whole organ. The QTM model is potentially applicable to other tracers and image modalities, such as DCE and DSC MRI, CT, and PET. The signal decay rate in Equation (3) should be correspondingly adjusted, that is, for contrast agent tracers that do not decay, and ln2/half‐life for radioligands used in PET.

To safely apply the proposed QTM, the assumptions in the development of it should be understood properly. First, there are several assumptions that should be understood carefully: The blood velocity is assumed to be the same as the velocity of the tracer. This is the basic assumption in the study of most tracer transport phenomenon. The pulsation of the blood velocity is ignored and is assumed to be constant in time because QTM is to measure an averaged effect of blood velocity and flow. The pulsation can be averaged out during a time interval that is longer than a cardiac period. In MRI, the ASL data were acquired within about 4 min; therefore, the data only reflect an average of blood velocity information. The assumption of impermeable vessel wall is reasonable because the tracers first transport in the renal vasculature before entering the renal tubules. The vessel wall was assumed to be impermeable, in other words, the model is basically one compartment model. For ASL data, the tracer signal decays at the rate of The tracer is on the arterial side and is generally regarded with little penetration through the vessel wall into the tissue compartment. This assumption serves as a reasonable start for this proof‐of‐concept work and for healthy subjects. Future work needs to incorporate permeability in the model for tissues that are peameable to tracers and for diseases that increase the vessel wall permeability. Two or more compartment model can be considered, and corresponding temporal convolution can be added to the transport equation of mass flux described in this work. The outlet flow was assumed to be equal for all the terminal capillary segments during the construction of the microvascular network. This is a reasonable assumption given that the tissue needs nutrition supply from the capillary and the capillary distribution in the renal cortex is quite homogeneous. This assumption can be removed in future work on the inhomogeneous case.

The proposed QTM has its limitations at this preliminary stage. First, microvasculature information is needed or at least the absolute blood volume value is required for accurate reconstruction of the blood flow map. The error in the estimation of blood volume may introduce errors in the final blood flow map. To produce a reasonable prediction of the absolute renal blood volume, one can use a population‐based value or patient‐specific predicted value. The population‐based renal blood volume may be obtained from large database and the patient‐specific renal blood volume may be calculated using MRI technique. Second, QTM uses ASL data with multiple PLDs, which means more data acquisition time is needed. The long data acquisition time for some patients may be intolerable. This problem may be solvable by shortening the data acquisition time for each PLD using Hadamard encoding. A navigator may be used to compensate motion artifacts from the long scan. , Third, QTM is a PDE based inverse problem in which the derivative of data plays a key role. Errors in the flow map may be caused by subject movement during data acquisition, which is an issue in abdominal imaging. However, motion artifacts can be reduced using state‐of‐the‐art motion correction algorithms and image registration algorithms. , Fourth, the microvasculature model of kidney employed here was from rat, and the computational power to calculate the forward problem in a large microvascular network limited the smoothness of the simulated data. Although there is no renal micrvasculature data for human available to our knowledge, the rat renal microvasculature is suggested to be very similar to that of human in both morphology and function. However, the renal vasculature in disease is different from that in health, and future studies on the renal microvasculature network of human in health and diseases are needed. Fifth, the simulation is based on transport equation and, therefore, favors QTM against Kety’s method, which is the so‐called inverse crime in the area of inverse problem. Fundamentally, Kety’s method assumes global arterial input for all voxels, ignoring the physics of local spatial variation contributing to local temporal variation, and any physics based modeling is biased against Kety’s method. In this work, the simulated data were first generated on a “microscopic” scale and then, voxelized and inverted on imaging voxel scale in QTM. This voxelization makes the resolution of QTM different from that of the forward problem, mitigating the inverse crime to some degree. Moreover, in vivo data provides a fair comparison between QTM and Kety’s method. Sixth, the non‐contrast MRA was used to obtain the prior information of main arteries of the kidney. We found that the renal artery structure for different subjects have a great similarity and it may be affected by the shape of kidney as well. We are thinking to train a convolutional neural network to predict the renal main artery structure given the shape of kidney to reduce the scanning time.

Several possible aspects of future work should be addressed to extend the perfusion quantification using QTM to clinical applications. We have demonstrated the feasibility of the proposed QTM in ASL data in the kidney. It may be applied to the other body regions including brain, liver, and breast, and the other types of tracer imaging methods including computed tomography (CT), PET, and single‐photon emission computed tomography (SPECT), as well as DCE MRI. Accurate quantification of regional blood flow, combined with quantitative susceptibility mapping, allows estimating the metabolic rate of oxygen consumption. Additionally, a flow phantom to experimentally validate the method with a known ground truth flow map would provide further validation. To get an accurate estimate of the flow map, QTM needs an accurate estimate of blood vasculature. Finally, the assumption of intravascular tracer distribution could be extended to the extravascular space in future work.

In summary, the proposed QTM is able to quantify blood perfusion by fitting the 4D time‐resolved image data to the voxelized transport equation with the aid of the reconstructed microvasculature information. QTM provides both flow and directional velocity maps and does not require an AIF.

FIGURE S1 (A) Illustration of microvasculature in one voxel, there are intersectional areas between vessels and voxel surfaces. The total intersectional areas on the six voxel surfaces are denoted as and , respectively. (B) illustration of voxelization and adjacent voxel relation. Negative x to positive x: Left (L) to Right (R), negative y to positive y: Posterior (P) to Anterior (A), negative z to positive z: Inferior (I) to Superior (S)

FIGURE S2 Illustration of boundary condition for the transport equation in cylinder

Click here for additional data file.