Tracer kinetic models as temporal constraints during brain tumor DCE‐MRI reconstruction

Purpose To apply tracer kinetic models as temporal constraints during reconstruction of under‐sampled brain tumor dynamic contrast enhanced (DCE) magnetic resonance imaging (MRI). Methods A library of concentration vs time profiles is simulated for a range of physiological kinetic parameters. The library is reduced to a dictionary of temporal bases, where each profile is approximated by a sparse linear combination of the bases. Image reconstruction is formulated as estimation of concentration profiles and sparse model coefficients with a fixed sparsity level. Simulations are performed to evaluate modeling error, and error statistics in kinetic parameter estimation in presence of noise. Retrospective under‐sampling experiments are performed on a brain tumor DCE digital reference object (DRO), and 12 brain tumor in‐vivo 3T datasets. The performances of the proposed under‐sampled reconstruction scheme and an existing compressed sensing‐based temporal finite‐difference (tFD) under‐sampled reconstruction were compared against the fully sampled inverse Fourier Transform‐based reconstruction. Results Simulations demonstrate that sparsity levels of 2 and 3 model the library profiles from the Patlak and extended Tofts‐Kety (ETK) models, respectively. Noise sensitivity analysis showed equivalent kinetic parameter estimation error statistics from noisy concentration profiles, and model approximated profiles. DRO‐based experiments showed good fidelity in recovery of kinetic maps from 20‐fold under‐sampled data. In‐vivo experiments demonstrated reduced bias and uncertainty in kinetic mapping with the proposed approach compared to tFD at under‐sampled reduction factors >= 20. Conclusions Tracer kinetic models can be applied as temporal constraints during brain tumor DCE‐MRI reconstruction. The proposed under‐sampled scheme resulted in model parameter estimates less biased with respect to conventional fully sampled DCE MRI reconstructions and parameter estimation. The approach is flexible, can use nonlinear kinetic models, and does not require tuning of regularization parameters.


INTRODUCTION
Dynamic contrast enhanced-magnetic resonance imaging (DCE-MRI) is a powerful technique that provides a quantitative measure of vessel permeability and interstitial volumes. In the brain, it characterizes the blood brain barrier (BBB) leakiness, which has proven to be valuable in several applications. 1 These include assessing conditions with large BBB breakdown such as gradation of brain tumors, 2,3 multiple sclerosis lesions, 4,5 and conditions with subtle and chronic BBB breakdown such as diabetes, 6 and Alzheimer's disease. 7 Outside the brain, DCE-MRI has applications in cancer assessment and therapeutic monitoring in several body parts including breast, 8,9 prostate, 10 and liver. 11 DCE-MRI involves a challenging trade-off between the achievable spatial resolution, temporal resolution, and volume coverage. Acceleration strategies that exploit redundancies along the time dimension have shown significant potential to improve these trade-offs. These include schemes such as view-sharing, [12][13][14] highly constrained back projection (HYPR), 15 and compressed sensing. [16][17][18][19][20][21] Several sparsifying spatio-temporal transforms have been proposed including spatio-temporal wavelet transform, spatio-temporal finite-difference, temporal Fourier transform. A major challenge with these "off-the-shelf" object models is that the modeling assumptions do not fit the data, which limits the achievable acceleration rates. Data-driven schemes that learn sparse representations from the data have been proposed, [22][23][24][25] and have shown to out perform off-the shelf transforms. However, these are often associated with highly nonconvex optimization. Furthermore, image reconstruction with existing transforms involves tuning one or more regularization parameters, which poses challenges to the standardization of these methods.
In this manuscript we explore the use of physical tracer kinetic models for constrained reconstruction. This approach has been used extensively in dynamic positron emission tomography (PET) imaging, [26][27][28][29] and has recently been adapted in MRI for the applications of relaxometry, [30][31][32][33] perfusion, 34,35 permeability, 36,37 and diffusion imaging. 38,39 Broadly, these methods can be classified into methods based on direct reconstruction of parameters from under-sampled data, or methods that use representations derived from parametric models as constraints in image reconstruction.
We propose a model-constrained approach for DCE-MRI reconstruction, where established contrast-agent kinetic models used in post-processing are employed as temporal constraints in reconstruction. From a specific kinetic model, and a physiological range of kinetic parameters, we construct a library of concentration vs. time profiles. Kinetic modelspecific temporal basis functions are derived from the library using the k-singular value decomposition (k-SVD) algorithm. 40 Through noiseless and noise-based simulations, we deduce a relation between the sparsity parameter in k-SVD and the complexity level of the kinetic model. We design a constrained reconstruction method where the kinetic modelbased temporal bases are used to constrain the recovery of concentration vs time profiles from under-sampled (k-t) data. We utilize an iterative multiscale optimization algorithm for improved robustness to undesirable local minima solutions.
The proposed approach has similarities with recent work on direct reconstruction of kinetic parameters from undersampled DCE-MRI data. 36,37 We use the same tracer kinetic model for reconstruction and post-processing to exploit the redundancy in the DCE-MRI pipeline. The major difference is in formulation of the optimization problem. Direct reconstruction involves estimation of kinetic parameters directly from under-sampled data. 36 When using the Patlak model, a Newton-based solver is used. When using the ETK model 37 a variable splitting strategy is used to iterate between sub-problems of data consistency, concentration time profile estimation, and kinetic parameter estimation. One major challenge is that the kinetic parameter estimation is treated as a black box. For models like ETK, this fitting problem is nonlinear and is applied to concentration vs. time profiles containing noise and potential artifact at every iteration. Errors in kinetic parameter estimation propagate to the main reconstruction step. A second limitation is substantial compute times as kinetic model estimation is performed at every iteration. In contrast, the proposed approach decouples kinetic parameter estimation from the reconstruction of concentration profiles. This avoids calling the computationally expensive kinetic parameter estimation during reconstruction. The proposed optimization iterates between data consistency, concentration profile estimation, and k-sparse projection of concentration profiles onto a set of temporal basis functions. Kinetic parameter estimation needs to be performed only once.
Since our formulation decouples reconstruction of concentration profiles from parameter estimation, it allows for flexibility to adapt to complex nonlinear kinetic models. Furthermore, since the sparsity parameter is fixed a priori, the proposed approach does not require any tuning of free parameters (e.g., regularization parameters). The flexibility allows for its potential utility in DCE-MRI of most organs and disease conditions. In this work, we demonstrate effectiveness with both the Patlak and extended Tofts-Kety (ETK) models, and demonstrate application to brain tumor assessment.

2.A. Tracer kinetic model-based temporal bases
A library of concentration vs. time profiles Z lxN is simulated using a kinetic model, an arterial input function (AIF), and a physiologic range of kinetic parameters (Fig. 1). l denotes the number of profiles in the library; and N denotes the number of time instances. For the ETK model, 41 we used the range: K trans ¼ 0 À 0:8min À1 in steps of 0:01min À1 , v p ¼ 0 À 60% in steps of 1%; v e ¼ 0 À 100% in steps of 1% to yield a library of size l Â N ¼ 494100 Â 50. Similarly, for the Patlak model, 42 we used the range: K trans ¼ 0 À 0:8min À1 in steps of 0:01min À1 , v p ¼ 0 À 60% in steps of 1% to yield a library size l Â N ¼ 4941 Â 50. We assume a hematocrit (hct) of 0.4, which is equivalent to the range of blood volume (v b ) between 0% and 100% as v b = v p / (1 -Hct). N was chosen as 50 to match our in-vivo DCE-MRI acquisition settings (i.e., temporal resolution of 5 s and the total scan time of 250 s). This can however be adjusted based on the temporal resolution and scan time of the DCE-MRI acquisition. A population-based AIF was used. 43 The settings of the Parker model that specifies the populationbased AIF were the same as described in Ref. [ 43 ]. The range of kinetic parameters was motivated by brain tumor DCE literature, 1 which suggests 0-0.34 min -1 for K trans , 0%-60% for v p assuming hematocrit of 0.4, and 0%-100% for v e . We expanded the K trans range by~2.5x, and used the full range for v p and v e to ensure conservative coverage of the kinetic parameter space. The k-SVD dictionary-learning algorithm 40 is then used to reduce the large library to a smaller dictionary of temporal basis functions (denoted by V rxN ). k-SVD represents any time profile in Z, for instance the pth row of Z, z p (t) as a sparse linear combination of basis functions v i t ð Þ from V: where r denotes the number of basis functions in V, and is chosen as r ¼ 100 ( l. q is the sparsity parameter. jju p jj 0 denotes the l 0 norm of the vector u p ¼ u 1 ; u 2 ; . . .; u r f g . z qÀsp p ðt) denotes the q-sparse projection of z p t ð Þ onto V. k-SVD jointly estimates the sparse coefficient matrix U lxr and the dictionary V rxN as: where u p denotes the pth row of U.

2.B. Image Reconstruction
We pose the estimation of the concentration vs time profiles C MxN (Mnumber of pixels; Nnumber of time frames), and the sparse coefficient matrix U Mxr from undersampled k-t space data (b) as: |fflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflffl ffl} data consistency ;s:t:;C ¼ UV;jju p jj 0 q;p ¼ 1;2;...;M f g ; |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} TK model constraint (3) C contains the concentration v.s time profile c x; t ð Þ for every pixel x 2 x; y ð Þ stacked row wise. c x; t ð Þ are constrained to be a q-sparse linear combination of the kinetic model-derived Þ denotes the forward model which maps C to the measured multicoil (k,t) data. F u denotes the Fourier Transform operator on a specified (k-t) under-sampling pattern. S m contains the receiver coil sensitivity maps. To estimate the coil maps, the standard sum of squares method is applied on high SNR multicoil images obtained after gridding reconstruction of the time collapsed raw k-t data; the coil maps are assumed to capture object phase. T is an operator that relates the concentration profile to the signal intensity profile s x; t ð Þ by the steady state spoiled gradient echo (SPGR) equation: where r 1 is the contrast agent relaxivity, TR is the repetition time, a is the flip angle, R 1 x; 0 ð Þ and M 0 x ð Þ are respectively the pre-contrast R 1 (reciprocal of T 1 ) and the equilibrium longitudinal magnetization. s x; 0 ð Þ is the pre-contrast first frame, which is fully sampled. The bracketed term in the second row of Eq. (5) resolves differences between the pre-contrast signal s x; 0 ð Þ and the predicted pre-contrast signal based on the baseline R 1 x; 0 ð Þ and M 0 x ð Þ maps (from a separate T 1 mapping acquisition). Similarly, the operation of mapping concentration profile from the signal intensity profile can be expressed as 44 : FIG. 1. Construction of dictionary of temporal basis functions from a specified tracer kinetic model (a). Based on a physiological range of kinetic parameters and an arterial input function (AIF), a library of concentration vs time profiles is generated (b). A subset of the profiles in the library is highlighted in red. Using k-SVD, the library is then reduced to a smaller set of temporal basis functions in a dictionary (c). The basis functions generated with the extended Tofts-Kety model is shown in (c). The basis functions themselves are not representative of kinetic model profile, and hence can be nonpositive. Instead, the linear combination of them is designed to mimic any profile in (b). Approximate MATLAB computational times respectively for generating the library (~400 000 profiles) and learning the dictionary were 11.5 min and 3.5 h.
[Color figure can be viewed at wileyonlinelibrary.com] Medical Physics, 47 (1), January 2020 We solve (4) by alternately (a) updating U using orthogonal matching pursuit (OMP) sparse projection, 40,45 and (b) updating C by enforcing consistency with acquired data. To be robust to spurious local minima, we use an iterative multiscale minimization approach, where we solve the problem at a coarser spatial resolution during the initial iterations and as the iterations proceed, we gradually update the resolution to its full resolution. This is achieved by multiplication of spatial Fourier Transform of s x; t ð Þ by a two-dimensional Gaussian filter (G k r ð Þ) specified by filter width k r ; where k r is initialized to 0.1% of k max , and gradually updated to 100 percent of k max , where k max specifies the extent of k-space coverage. This heuristic strategy is used in several nonconvex problems such as in image registration, 46 and recently in MR-fingerprinting. 47,48 Starting with an initial guess obtained from ð Þ, we iterate until a stopping criterion of

2.C. Simulations
The sparsity parameter q in Eq. (2) is determined based on simulation studies with the Patlak and the ETK models. Noiseless simulations are performed and the mean approximation error l err ¼ 1 l P l p¼1 jjz p t ð Þ À z qÀsp p t ð Þjj 2 2 , and maximum approximation error: Noise-based simulations were performed for broad ranges of kinetic parameter values to (a) determine any systematic bias and uncertainty in the kinetic parameter space that may be induced by sparsity-based modeling of the concentration time profiles, and (b) to deduce the correspondence between the sparsity level (q) and the kinetic model.
Noisy concentration profiles were obtained as: where n t ð Þ denotes i.i.d. white Gaussian noise with zero mean and 0.005 standard deviation, which was chosen to match the typical signal-to-noise ratio (SNR) from in vivo brain DCE-MRI data acquired at our institution on a 3T commercial system with an eight-channel head array coil. Noise in the concentration time profiles was assumed to be additive i.i.d. Gaussian. Concentration time profiles are real valued (negative values can occur in the presence of noise), and have a one-to-one mapping with real-valued signal intensity time profiles in the forward model. Monte-Carlo simulations with 500 realizations of n t ð Þ were performed to evaluate the bias and uncertainty in estimating kinetic parameters from (a) the noisy profiles z n p t ð Þ, and (b) the qsparse projections of z n p t ð Þ on V: z n;qÀsp p t ð Þ. We performed covariant error analysis for two parameters (K trans , v p ) with the Patlak and the ETK model over a broad range of kinetic parameters. With both the models, we evaluated the bias and uncertainty in estimating K trans and v p before and after q-sparse projections. With the ETK model, for simplicity, we focus only on analysis in a two dimensional space with a fixed v e = 0.6. The open-source Rocketship package 49 was used for kinetic parameter estimation.

2.D. Evaluation with a digital reference object
An anatomically realistic brain tumor DCE-MRI digital reference object (DRO) was generated based on the method and data described in Ref. [ 50 ]. Briefly, the population-based AIF with the Parker model, known kinetic parameters, the ETK model, and the steady state spoiled gradient signal equation was used to generate the dynamic images. We then multiplied by coil sensitivities, took the Fourier Transform, and added realistic complex Gaussian noise to each channel. Coil maps, noise covariance matrix, and the signal to noise (SNR) level were obtained from in-vivo data acquired at 3T. Comparisons were performed at a SNR = 30 to mimic measurements at 3T.
This phantom data was retrospectively under-sampled using a randomized golden-angle Cartesian (GOCART) sampling pattern, 51 and evaluations in fidelity of the kinetic parameters were performed at under-sampling factor of R = 20. GOCART 51 is originally a three-dimensional (3D) golden angle Cartesian sampling scheme, with random sampling of the ky-kz phase encode locations along each Cartesian radial spoke. In this study, we perform retrospective under-sampling in the kx-ky plane in a representative slice. This strategy was chosen to simulate k y -k z under-sampling in prospective acquisitions.
Reconstruction was performed using the fully sampled direct inverse Fourier Transform-based approach (considered as reference), and the proposed constrained reconstruction approach with full-sampling (R = 1), and under-sampling (R = 20). Kinetic modeling was performed both using the ETK model, and a simpler Patlak model. The latter was used to analyze any errors due to model mismatch. The proposed constrained reconstruction was implemented with temporal dictionaries derived from the kinetic model (Patlak or ETK) that corresponded to the same kinetic model used for subsequent parameter estimation.

2.E. Evaluation with in-vivo data
We reviewed 110 fully sampled DCE-MRI raw datasets from patients with known or suspected brain tumor, receiving a routine brain MRI with contrast on a clinical 3T scanner (HDxt, GE Healthcare, Waukesha, WI). The data acquisition was based on a 3D Cartesian SPGR sequence with field of view (FOV): 22 9 22 9 4.2 cm 3 ; spatial resolution: 0.9 9 1.3 9 7.0 mm 3 ; temporal resolution: 5 s; 50 time frames; and eight receiver coils; flip angle of 15°, TE/ TR = 1.3 ms/6 ms. Driven equilibrium Single Pulse Observation of T1 (DESPOT1) was performed before the DCE sequence, where three images with flip angles of 2°, 5°, 10°w ere acquired to estimate T1 and M0 maps before the contrast arrival. Gadobenate dimeglumine (Multihance, Bracco) (0.05 mmol/kg) was administered into an upper extremity vein using a power injector (ACIST EmpowerMR Injector, Bracco), at a rate of 3 ml/s, followed by a 20 ml saline flush.
Of these 110 cases, we identified a cohort of 12 cases, which had different brain tumor characteristics (shape, size, heterogeneity), and also had enhancing tumors of atleast 1 cm (as determined by standard bi-directional assessment). 52 The demographics of these patients are shown in Table I, and the post-contrast images (last spatial frame from the DCEscans) are shown in Fig. 2. The protocol was approved by our institutional review board (IRB).
(k-t) under-sampling was performed retrospectively on fully sampled raw data using the GOCART (randomized golden angle Cartesian) sampling trajectory 51 at acceleration factors R = 20 and R = 40. Image reconstruction was performed with the proposed dictionary based approach, an existing compressed sensing approach that uses a temporal finite difference (tFD) sparsity constraint, 19 and compared against the reference fully sampled inverse Fourier Transform based approach. ETK derived bases with a fixed sparsity level of q = 3 was used in the proposed approach. The ETK model was chosen as it accounts for backflux of contrast from the extravascular space to the plasma, which in turn improves the accuracy of K trans estimation, and has shown to be applicable to brain tumor data. 53 All the patient datasets were acquired with fixed injection timing, however timing delays between 5 and 10 s (1-2 frames) existed amongst different patients.  2. Post-contrast images of 12 brain tumor cases with different brain tumor characteristics (shape, size, heterogeneity). All cases had enhancing tumors of at least 1 cm as determined by standard bi-directional assessment. 52 Fully sampled raw multicoil (k-t) space data from these patients were used as reference in retrospective under sampling studies. Tofts-Kety parameter estimation was performed in the tumor regions of interests as marked by the red shaded regions. [Color figure can be viewed at wileyonlinelibrary.com] As described earlier, a population based AIF with a fixed delay was used to generate the library. Patient-specific AIF delays were estimated as described by Lebel et al. 54 Briefly, the k-space origin was frequently sampled, and plotted as a function of time. The region of maximum slope was regressed to the baseline to determine the bolus arrival time.
Either padding zeroes initially to the acquired data or omitting the last time frames corrected for any delay mismatch to the library. The tFD-based formulation is convex and is guaranteed to achieve the global minimum. Therefore the multiscale optimization heuristic was not applied during tFD optimization. We used the alternating direction method of multipliers (ADMM) algorithm where the stopping criterion was if the rate of change between reconstructions at successive iterations fell below 10 À5 percent. The regularization parameter in tFD constrained reconstruction was tuned to provide the smallest normalized root mean squared image reconstruction error (nRMSE) in tumor ROIs with respect to the reference fully sampled datasets. All reconstructions were implemented in MATLAB (The Math-Works, Inc., Natick, MA) and executed on an Intel core i7 3.5 GHz machine with 32 GB memory. The convergence of the proposed multiscale iterative optimization was evaluated empirically. Reconstruction estimates with different initializations of the concentration profiles were compared: zero filled reconstruction (C init ¼ A H b ð ÞÞ; low spatial resolution estimate obtained from the center 3x3 window of the k-space data in every time frame (C init ¼ C low:res Þ; and from the reference fully sampled data After image reconstruction, the ETK model was used to estimate the kinetic parameters with a population based AIF. 43 Bland-Altman analysis was performed to evaluate systematic bias and uncertainty of the reconstructed b K trans andv p maps (from the proposed and tFD approaches) with respect to the reference fully sampled K trans ref and v p;ref . Comparisons using v e maps were not considered, as its estimation is associated with high uncertainty with the ETK model. 49 To evaluate error in the kinetic maps on using the ETK model to constrain the reconstruction of time intensity curves, we analyzed the kinetic maps from the proposed reconstruction against conventional direct inverse Fourier Transform reconstruction on fully sampled data (R = 1).
On one of the in-vivo datasets where the Patlak model produced less than one percent modeling error with the population based AIF, the proposed reconstruction scheme implemented with the Patlak dictionary was also compared against the direct Patlak parameter estimation method. 36 The direct estimation approach was implemented based on open source code (https://github.com/usc-mrel/DCE_direct_ recon). The reconstruction quality in the resulting K trans and v p maps from under-sampled data (at R = 10-30) were compared against the maps obtained by conventional fully sampled (R = 1) inverse Fourier Transform reconstructions. Figure 3 shows the maximum and average approximation errors (max err ; and l err )between the concentration vs time profiles in the library, and the profiles obtained from q-sparse projections onto V at different sparsity levels (q). q-sparse projections of the curves generated from the Patlak and the ETK models are respectively shown in Figs. 3(a), and 3(b). These curves are chosen to represent different types of tumor enhancement dynamics. 55 With q = 1, we observe considerable bias in approximating the kinetic model generated curves with both the Patlak and the ETK models. However, for the Patlak model, a choice of q ≥ 2 provided excellent agreement with the profiles in the library (max err =l err ¼ 10 À28 %=10 À30 %Þ. Similarly, for the ETK model, a choice of q ≥ 3 approximated the profiles in the library with (max err =l err ¼ 2%=0:008%Þ. Figure 4 and 5 demonstrates the bias and uncertainty in estimating kinetic parameters in presence of noise. Over a broad range of kinetic parameters, we observe that estimating the kinetic parameters from noisy profiles and the q-sparse projected profiles are equivalent when q ≥ 2 (for the Patlak model), and q ≥ 3 (for the ETK model). Based on these simulations, we fixed q = 2 for the Patlak model, and q = 3 for the ETK model. Figure 6 shows evaluations on the brain tumor DRO which is constructed based on the ETK model. When the simpler Patlak model is applied on the reconstructed concentration time profiles, an under-estimation in K trans (a factor of about 10 fold) and over estimation in v p (a factor of about 2.5 fold) is observed. This bias is observed in both the reference fully sampled inverse Fourier Transform reconstruction, and the proposed reconstruction suggesting that model selection error is independent of the choice of the reconstruction [see Fig. 6(b)]. When evaluated against the kinetic parameter estimates from reference reconstructions, the proposed reconstruction showed robust quality maps from 20 fold under-sampled data [Figs. 6(a) and 6(b)]. This is also depicted in the difference maps where the reconstruction error lies at the level of background noise. Figure 7 shows the kinetic maps from the proposed reconstruction against conventional direct inverse Fourier Transform reconstruction on fully sampled data (R = 1). Two representative brain tumor datasets that have different spatial characteristics are shown: (a) glioblastoma with thin rim, necrotic core, connected to a solid tumor; (b) meningioma with homogenous spatial tumor characteristics. As depicted in the kinetic maps, the maps obtained from the two approaches are qualitatively equivalent, with the error being in the background noisy regions. This is also highlighted in the difference maps. For both the reconstructions, the ve maps are noisy due to the increased uncertainity in estimating ve from short scan times (5 min in this study). 56 Figure 8 shows the evolution of the objective function in Eq. (3) as a function of CPU reconstruction time with different initializations of concentration time profiles (a) from low-resolution dynamic images C low:res: ; (b) from zero-filled dynamic images C ¼ A H b ð Þ; (c) from reference dynamic images C ref . The multiscale optimization gradually updates the complexity of the problem. Due to spatial low-pass filtering, the under-sampling artifacts in initial iterations are considerably reduced making the problem well-posed. b C is updated gradually with increasing resolution, as a result of which a monotonic convergence is observed. We empirically found this approach to be robust to local minima; the final solutions were identical with different initializations. Figure 9 shows retrospective under-sampling comparisons of K trans at R = 20. tFD reconstructions resulted in considerable under estimation of K trans in 9 of the 12 cases, while the proposed method was found to be robust to this bias. tFD also relied on adjusting the regularization parameter. In contrast, the proposed parameter free reconstruction provided K trans estimates closer to that of the reference. It also provided superior fidelity in maintaining spatial characteristics of the tumors in all cases (e.g., depiction of thin tumor boundaries in cases 1 to 5). Figure 10 shows Bland-Altman plots of the difference between estimated TK parameters (at R = 20, R = 40) and . The proposed reconstruction employs the ETK-based dictionary in (i), and the Patlak-based dictionary in (ii). From (i) and (ii), it is seen that when a simpler Patlak model is used in post-processing, the resulting K trans is under estimated in comparison to the ETK-based K trans estimates (about a factor of 10-fold). This under estimation is observed in all the reconstructions (a-c) suggesting kinetic model selection error is independent of the type of the reconstruction scheme. From (a-c), it can be seen that proposed reconstruction shows good quality in the kinetic parameter estimates at R = 1 and R = 20, which is also highlighted in the difference maps scaled by a factor of 3 in (d-e). [Color figure can be viewed at wileyonlinelibrary.com] FIG. 7. Comparison of kinetic parameters obtained from concentration time profiles from reference fully sampled data (first row) against the profiles after 3sparse projections onto the extended Tofts-Kety (ETK) dictionary (second row). The third row denotes the difference map that highlights the kinetic modeling errors. The difference map is scaled up by a factor of 3 for better visualization. Two representative brain tumor datasets that have different spatial characteristics are shown: (a) glioblastoma with thin rim, necrotic core, connected to a solid tumor; (b) meningioma with homogenous spatial tumor characteristics. The kinetic maps in the first two rows are observed to be qualitatively equivalent suggesting 3-sparse projection mimics ETK modeling. Compared to the K trans , and v p maps, the uncertainty in v e is large attributed to the short scan time duration of 5 min. [Color figure can be viewed at wileyonlinelibrary.com] reference TK parameters on all the 12 cases combined. In comparison to tFD, the proposed approach showed reduced bias and reduced certainty during estimation of K trans , and v p . A systematic bias of under estimating K trans , and v p was present in tFD, which is also qualitatively shown in Fig. 8. Figure 11 shows the comparison of the estimated TK parameters from the proposed approach with Patlak dictionary against the direct Patlak parameter estimation at R = 20, 30, 40 on case 7. Both the proposed approach and the direct reconstruction provided good spatial fidelity of the TK maps as depicted in Fig. 11(i). The proposed approach however demonstrated subtle benefits in terms of reduced bias, and reduced noise amplification as shown in the Bland-Altman plots of Fig. 11(ii).

DISCUSSION
We have developed a new DCE-MRI reconstruction approach that applies kinetic models routinely used in postprocessing as temporal constraints during reconstruction. Based on simulation studies, we deduced a relation between sparsity parameter q in k-SVD to the complexity of the kinetic model. We have demonstrated equivalence of Patlak and ETK models with dictionaries constructed respectively with q = 2 and q = 3. This approach exploits the smooth time intensity DCE patterns by using temporal basis functions derived from a kinetic model. This is in contrast to generic off-the-shelf transform bases that are blind to the kinetic model behavior of the time intensity profiles. We also proposed a robust . tFD also relied on tuning of a regularization parameter. In contrast, the proposed parameter-free model-based reconstruction provided K trans estimates closer to that of the reference, and has improved fidelity in preserving spatial characteristics of the tumors (e.g., thin boundaries of the tumor, see arrows in cases 1-5). [Color figure can be viewed at wileyonlinelibrary.com] multiscale iterative optimization algorithm to solve the resulting l0 norm based nonconvex objective function. We empirically demonstrated robustness to local minima. The tFD approach has a convex formulation with a guaranteed global minimum, and therefore, no multiscale optimization heuristics were applied during tFD minimization. In-vivo validation with 12 brain tumor cases demonstrated superior recovery performance with the proposed method compared to tFD (reduced bias, uncertainty in kinetic mapping, and better spatial fidelity of kinetic maps) at up to R = 40.
The proposed framework can be extended in several ways. A uniform grid of kinetic parameters was used in this study to generate the library of possible concentration profiles from a chosen kinetic model. However, it is possible to perform application-specific discretization of the kinetic parameters to improve sensitivity and accuracy in modeling time curves that lie in a particular zone in the kinetic parameter space. We have demonstrated that for the 2-parameter Patlak model, a choice of two temporal bases from the dictionary was adequate to reliably model the concentration temporal profiles. For the 3-parameter ETK model, three temporal bases were adequate. We expect that if this approach is extended to more complex models (e.g., fast exchange, shutter speed, two compartment exchange model), the complexity of the bases representation may also increase. Increasing model complexity is expected to place more stringent limits on the maximum achievable acceleration rates.
This study demonstrates the utility of the kinetic modelbased reconstruction approach in the application of brain tumor imaging. The framework may be extended to DCE-MRI of other diseases and body parts by appropriately considering different ranges of kinetic parameters and kinetic models. Other applications would also warrant validation with digital reference objects specific to that body part, disease, and kinetic parameter ranges. Complementary constraints such as spatial sparsity could be added to further improve the recovery.

4.A. Limitations of the study
A limitation of this feasibility study is the use of population-averaged AIF. 43 Population-averaged AIFs are known to produce a potential bias in the final kinetic maps, 57 however in this study, the bias identically affects the reference maps, and maps produced by the proposed reconstruction and the temporal finite difference reconstruction that is used for comparison. The framework can be extended to account for patient-specific AIFs. For instance, the full patient-specific AIFs can be obtained from a preprocessing reconstruction or could be characterized by including richer dictionaries that parameterize the shape and amplitude of the AIF. 58,59 Our preliminary findings show that with whole brain scans and transform sparsity reconstruction, we could extract good fidelity AIFs that are free of inflow enhancement artifacts. 54 Such a reconstruction could potentially be a preliminary pre-processing step. These extensions are a scope of our future work.
In this study, our DCE-MRI protocol employed a flip angle of 15°, and in-vivo data were acquired at ½ dose (0.05 mMol/kg). With these settings, there exists some nonlinearity in mapping between concentration profiles and signal intensity time profiles at concentrations > 1.0 mM. The linearity can however be improved with the use of flip angles ≥ 25°. R2* effects were not included in the forward imaging model. 60 Our DCE-MRI scans were performed with ½ dose (0.05 mmol/kg), and used a short TE of 2 ms. We have examined several clinical datasets at our institution and have found phase and R2* effects to be insignificant in tissue and in vessels. FIG. 10. Bland-Altman plots of (a) the difference between estimated K trans (at R = 20, R = 40) and reference K trans ; (b) the difference between estimated vp (at R = 20, R = 40) and reference vp; for the proposed (left column) and tFD (right column) reconstructions. Each dot corresponds to one pixel within the tumor ROIs of all the 12 cases. The mean and 1.96 times the standard deviation (l AE 1.96r) of the difference entities are quantitatively shown. These are also qualitatively marked by the solid red and dotted red lines. As seen from the plots, the proposed approach had lower bias (l) and uncertainty (r) in estimating K trans , and v p in comparison to tFD. tFD depicted a systematic bias in underestimating K trans , and v p in comparison to the proposed approach This can also be noted from the qualitative comparisons in Fig. 9. [Color figure can be viewed at wileyonlinelibrary.com] We therefore did not consider R2* or off-resonance effects, but these could be easily added to the forward model. v e has shown clinical value in studies with long scan times of 10 minutes or more (e.g., Ref. [ 61 ]). At our institution (and most of our peer institutions) the standard of care brain tumor protocol utilizes a DCE-MRI scan time of 5 min, which is insufficient to recover v e . This is documented in the literature where estimation of v e had high uncertainty with short scan times, 56 and longer scan times o upto 10 min are recommended. 56,62 We have included a range of values for v e in the library because accounting for backflux is known to improve the estimation of K trans and v p . 41 In this study, we have only focused our analysis only on tumor ROIs. We intend to perform comprehensive statistical studies against controls (nonleaking areas from the whole FIG. 11. (i) Comparison of K trans and v p derived from the methods of direct reconstruction of Patlak parameters 36 , and the proposed approach with a Patlak dictionary at R = 10, 20, 30 against the reference fully sampled reconstruction followed by Patlak modeling (R = 1). (ii) shows the Bland-Altman plots of the difference between estimated K trans (at R = 10, 20, 30) and the reference K trans for the direct reconstruction and proposed reconstruction with the Patlak dictionary. Each dot corresponds to one pixel within the tumor ROI of case 7 as marked in Fig. 2. The mean and 1.96 times the standard deviation (l AE 1.96r) of the difference entities are quantitatively shown in (ii). Both direct reconstruction, and the proposed dictionary reconstruction provides good spatial fidelity in the reconstructed maps (in i). From the Bland-Altman plots in (ii), the proposed approach demonstrated subtle improvements over the direct reconstruction in terms of a smaller bias and reduced noise amplification. [Color figure can be viewed at wileyonlinelibra ry.com] brain) in a future study with prospective under-sampling that can enable whole brain coverage.
Comparison of the proposed approach against direct reconstruction of Patlak parameters showed subtle improvements in estimating the TK parameters with less noise amplification, and reduced bias on a single dataset. Comprehensive evaluation against the direct reconstruction method when using nonlinear models and multiple datasets were not performed in this study as this warrants establishing a detailed study of optimization routines in direct reconstruction (e.g., Newton based, 36 variable splitting 37 ). An important distinction of the proposed work is that it decouples computationally expensive kinetic parameter estimation from the reconstruction of concentration vs time profiles. Kinetic parameter estimation need only be performed once as a final step. In contrast, direct reconstruction, as described in Guo et al., 37 requires kinetic parameter estimation at every iteration. This itself can require on the order of minutes to hours for typical DCE-MRI datasets with sizes of at least 128 9 128 pixels and 50 time frames. Precise reconstruction times depend on the choice of optimization solver and their implementation (e.g., efficiency and whether they exploit GPUs). We have not performed detailed comparisons of reconstruction times, as this is beyond the scope of this feasibility study.
The T1 maps were estimated prior to reconstruction using DESPOT1 with three flip angles. 63 However, using fully sampled data, the joint estimation of T1 and kinetic parameter maps has recently been shown to improve accuracy of DCE kinetic parameter maps. 64 An extension of the proposed framework to include joint T1 estimation would warrant the inclusion of multiple T1-based simulated concentration curves in the library, and exploration of superior learning approaches (alternate to k-SVD) that offer better compression capabilities to efficiently represent a richer library. Furthermore, there could be a number of approaches to improve pre-contrast T1 mapping in a separate step. For instance, increasing the number of flip angle measurements, use of constrained imaging methods (e.g., model-based reconstruction, MR fingerprinting).
The proposed approach has similarities and important distinctions with prior art. Similar to MR-fingerprinting, 48,65 our approach exploits physical models for reconstruction. However, it does not modify the acquisition parameter settings. It takes a two-step approach of first reconstructing the concentration time profiles, and then estimating the kinetic parameters in a final step. In comparison to MR-Fingerprinting, our approach is sensitive to motion because the basis functions do not account for motion. However, if reasonable estimates of the motion deformation fields are known or can be estimated from the data, it can be corrected by integration into the forward model. 66,67 Our experiments in this work show incoherence along time benefits the reconstruction. However, a detailed evaluation of various sampling pattern choices including coherent sampling, and evaluation against incoherent sampling is yet to be done, and is a scope of future work.
Data inconsistencies such as motion, or B1 non-uniformity, may violate the assumption of the appropriateness of the kinetic model on the concentration time profiles. This is equally true with existing compressed sensing methods. However, the framework can seamlessly accommodate prior information in the forward model to improve data consistency (e.g., integration of motion maps, B1 maps). The proposed reconstruction assumes the chosen kinetic model to be appropriate to the data. While the kinetic model of choice can be motivated based on the application at hand, the framework is flexible to generate comprehensive libraries from more than one kinetic model. Future application-specific studies with chosen kinetic models are however needed to deduce the relation between the complexity of the library and the sparsity parameter (q) during dictionary generation, and subsequently acceleration capabilities.