A Semi-parametric Technique for the Quantitative Analysis of Dynamic Contrast-enhanced MR Images Based on Bayesian P-splines

Dynamic Contrast-enhanced Magnetic Resonance Imaging (DCE-MRI) is an important tool for detecting subtle kinetic changes in cancerous tissue. Quantitative analysis of DCE-MRI typically involves the convolution of an arterial input function (AIF) with a nonlinear pharmacokinetic model of the contrast agent concentration. Parameters of the kinetic model are biologically meaningful, but the optimization of the non-linear model has significant computational issues. In practice, convergence of the optimization algorithm is not guaranteed and the accuracy of the model fitting may be compromised. To overcome this problems, this paper proposes a semi-parametric penalized spline smoothing approach, with which the AIF is convolved with a set of B-splines to produce a design matrix using locally adaptive smoothing parameters based on Bayesian penalized spline models (P-splines). It has been shown that kinetic parameter estimation can be obtained from the resulting deconvolved response function, which also includes the onset of contrast enhancement. Detailed validation of the method, both with simulated and in vivo data, is provided.

. For DCE-MRI, after the administration of a contrast agent, such as Gadolinium diethylenetriaminepentaacetic acid (Gd-DTPA), a dynamic imaging series is acquired. Typically, T 1 -weighted sequences are used to assess the reduction in T 1 relaxation time caused by the contrast agent. The contrast agent concentration time series can be estimated from the observed signal intensity using proton density weighted MRI after calibration or using multiple flip angle sequences [3], [4]. With Gd-DTPA, the agent does not enter into cells, so DCE-MRI depicts exchange between the vascular space and the extra-vascular extra-cellular space (EES).
Current approaches to quantitative analysis of DCE-MRI generally rely on non-linear pharmacokinetic models.
These models are usually derived from the solution to a system of linear differential equations, which describe the blood flow in the tissue [5]. In practice, single-compartment models may not always be suitable, and thus more complex models have been suggested. They include the "extended" Tofts-Kermode model (see Eqn. 2) and the tissue homogeneity approach [6], [7]. However, non-linear regression models are difficult to optimize and the estimation of the parameters depends on the initial values of the algorithm [6].
Recently, model-free techniques have received extensive attention in quantitative imaging. Neural networks are used for tissue classification [8], [9] and semi-parametric methods in the kinetic modeling of dynamic PET imaging [10], [11]. For the quantification of first-pass myocardial perfusion, Jerosch-Herold et al. [12] proposed a modelfree approach. The formulation of the response function based on a B-spline polynomial representation, however, is ill-conditioned and a first-order difference penalty spline (Tikhonov regularization) has been imposed.
As an alternative, Bayesian inference has also been investigated to replace traditional least-square fitting algorithms in MRI [13], [14]. Bayesian methods allow a more accurate description of the estimation uncertainty and can effectively reduce bias [15]. Furthermore, Bayesian hierarchical models can incorporate contextual information in order to reduce estimation errors [16]. A popular, general approach for semi-parametric modeling is based on Penalty splines, or P-splines [17], [18]. The function under consideration is approximated by a linear combination of a relatively large number of B-spline basis functions. A penalty based on k-th order differences of the parameter vector ensures smoothness of the function. For selecting the penalty weight (or smoothing parameter), the L-curve method [19] or cross validation can be used. In Bayesian frameworks, P-splines regression parameters can be estimated jointly with the penalty weight and hence allow for adaptive smoothing [20], [21].
In this paper, we propose a Bayesian P-spline model to fit the observed contrast agent concentration time curves.
The model uses a locally adaptive smoothing approach as the observed time series signal varies rapidly in the first minute after injecting the contrast bolus. The proposed algorithm provides a semi-parametric deconvolution approach and results in a smooth response function, along with the corresponding estimates of uncertainty. Kinetic parameters are robustly derived by fitting a non-linear model to the estimated response function. In addition, the exact onset of the contrast uptake can be determined by using information derived from the Bayesian estimation process. Detailed validation of the method both with simulated and in vivo data of patients with breast tumors is provided.

A. Standard kinetic models for DCE-MRI
The standard parametric model for analyzing contrast agent concentration time curves (CTCs) in DCE-MRI is a single-compartmental model [22], where the solution can be expressed as a convolution of the arterial input function (AIF) with a single exponential function [23], i.e., Here, C t (·) denotes the concentration of the contrast agent, C p (·) denotes the AIF, K trans represents the volume transfer constant between blood plasma and EES, and k ep represents the rate constant between EES and blood plasma. In this study, we use an extended version of the Tofts-Kermode model [6] as a reference parametric model given by where the additional parameter v p represents the fraction of contrast agent in the vascular compartment of the tissue.
The arterial input function C p describes the input of the contrast agent into the tissue. DCE-MRI studies typically use a standardized double exponential AIF given by [24] with values a 1 = 24 kg/l, a 2 = 6.2 kg/l, m 1 = 3.00 min −1 and m 2 = 0.016 min −1 [25]. In Eqn. 3, D is the actual dosage of tracer in mmol/kg. With this explicit form of the AIF, the convolution in Eqn. 2 can be derived analytically, i.e., 1) Discrete Deconvolution: Mathematically, a more general expression of Eqn. 1 can be written as: where f (t) is the response function in the tissue. Assuming that C p (·) and f (·) are constant over small intervals ∆t, a discretized form of Eqn. 5 is given by where C t is measured on discrete time points τ 1 . . . , τ n . The n×T matrix A may be interpreted as a convolution operator and is defined by where n i is the maximum index j for which τ i ≤ t j holds. It is worth noting that the input function C p (t) is measured -or evaluated from Eqn. 3 -at time points t 1 , . . . , t T , which can be different from τ 1 , . . . , τ n .
2) Penalty Splines: By solving Eqn. 6, the response function f (t) can be deconvolved from C p (t). However, this system may be numerically unstable, i.e., the deconvolved response function is susceptible to noise. To overcome this problem, we assume that f (t) is a smooth, k-times differentiable function. To this end, a B-spline representation of the response function is used in this study, i.e., where B is the n×p design matrix of kth order B-splines with knots s 1 , . . . , s p+k . In vector notation f = (f (t 1 ), . . . , f (t T )) and Eqn. 8 may be expressed as Accordingly, Eqn. 6 can be written as where D = AB is a n×p design matrix, representing the (discrete) convolution of the AIF with the B-spline polynomials.
To enhance the numerical stability, a penalty on the B-spline regression parameters is introduced, such that These models are known as penalty splines or P-Splines [17], [18] as they penalize the roughness of the function f (t), and therefore act as a denoising method. As C t exhibits a sharp initial increase followed by a sharp decrease at the beginning of the dynamic series, the penalty has to be locally adaptive. To this end, we use a Bayesian hierarchical framework for parameter inference.
3) Bayesian hierarchical framework: In Bayesian inference, a priori information, i.e., information available before observation of measurement, has to be expressed in terms of probability distributions. Here, our prior knowledge is the assumption about the model. In this paper, we assume that the observed contrast agent concentration C t is noisy realization of the true model (Eqn. 10), i.e., where D t denotes the tth row of D. That is, a priori the error is assumed to be Gaussian distributed with unknown variance σ 2 . We use a relatively flat prior for the variance parameter, i.e., where IG denotes the Inverse Gamma distribution.
The penalty in Eqn. 11 can be expressed as a priori distribution on β [21], where δ 2 t is the variance of e t . For a locally adaptive estimation of δ 2 t , the following prior model is used which allows time varying smoothness penalties.

4) Evaluation of the posterior distribution:
Bayes' theorem was used to calculate the posterior distribution of the parameter vector which is, up to a normalizing constant, given by A closed-form solution of Eqn. 16 is not possible, and thus, Markov chain Monte Carlo (MCMC) techniques have been used to assess the posterior distribution [29]. The full conditional of β, i.e., the joint distribution of β given all other parameters and the data, C t , is a p-dimensional multivariate normal distribution where R is the inverse covariance matrix of the prior distribution of β [21]. The full conditionals of both variance parameters are independent Inverse Gamma distributions To assess the posterior distribution, samples of the parameters β, σ 2 and δ 2 are drawn alternately from Eqn. 17, Eqn. 18 and Eqn. 19. After a sufficient burn-in period, the MCMC algorithm produces samples from the posterior distribution.
C. Estimating kinetic parameters from the semi-parametric technique The advantage of parametric models defined by Eqns. 1 and 2 is that they contain biologically meaningful and interpretable parameters. The semi-parametric technique only provides a fit to the CTC and a de-convolved and de-noised response function, but not kinetic parameters. In this section, we introduce methods for estimating kinetic parameters and the onset of the contrast uptake from the estimated response function. We make use of the fact that the Bayesian approach provides a rich source of information via the posterior distribution of the response function.
1) Determining the onset of contrast uptake: In practice, the time delay between the injection of the contrast agent and the arrival of the tracer in the tissue of interest is unknown. However, it is important to correctly estimate the delay for a robust estimation of the kinetic parameters [30]. The Bayesian approach yields information on the uncertainty for each parameter in the model and for all transformations of the parameters such as the estimated contrast agent concentration C t = Dβ. This results in a point wise Credible Interval (CI) for the contrast agent concentration. From this, we can compute the minimal time t * , where the 99% CI does not cover 0. That is, at time t * we are 99% confident that the contrast agent has already entered the tissue.
Assuming that the initial slope of the contrast agent concentration time curve is approximately linear, we can compute the onset time by drawing a line from C t (t * ) to 0 with the gradient dC t (t * )/dt. The first order derivative Since f (t) is a spline, the derivative can be computed as [31] where B (k−1) is the design matrix of (k − 1)th order B-Splines and γ j is defined via In this paper, we propose the following algorithm to estimate the onset of the contrast uptake: 1) Find the minimum time t * , for which the contrast concentration significantly exceeds zero, 2) Compute the gradient of the estimated CTC at t * , 3) Calculate the enhancement onset time as DCE-MRI studies typically assume that the onset of the enhancement is the same over the whole region of interest (ROI). In this case, the median of the estimated t 0 for all voxels in the ROI may be used as an estimate of the global value. However, for larger ROIs, local estimates of the onset may be required, and can be computed from the proposed technique.
2) Obtaining kinetic parameters: The semi-parametric technique provides a de-convolved and de-noised response function. In order to obtain kinetic parameters, we can fit a non-linear function to the estimated response. To this end, the following model has been used: where T c = v p /K trans is the transit time through the capillary, v e = K trans /k ep is the volume fraction of EES, E is the extraction fraction, and F p is the mean plasma flow. This model is similar to the adiabatic approximation of tissue homogeneity (AATH) [7]. In this model extraction fraction E and mean plasma flow F p may be not identifiable, but the product EF p = K trans is.
It should be noted that the Bayesian methodology provides not just one response function, but a (posterior) distribution of response functions. In order to obtain a distribution of the estimated parameters, the model in Eqn. 20 can be fitted to each response function in the sample using the Levenberg-Marquardt optimization algorithm. The median of the posterior distributions may be used to estimate the parameters. The estimation error can be computed using the standard errors across the samples and intervals can be constructed from quantiles of the posterior distributions.

D. Data Collection
Simulated DCE-MRI data of CTCs with known kinetic parameters were previously published in [6]. We use the first series which was designed to be representative of data acquired from a breast tumor. Data was simulated using In vivo data was derived from twelve patients with primary breast cancer (median age 46 years; range 29-70).
Each patient was scanned twice, once before and once after two cycles of chemotherapy. Scans were performed with a 1.5 T Siemens MAGNETOM Symphony scanner (TR = 11 ms and TE = 4.7 ms; 40 scans with four sequential slices were acquired in about 8 minutes). A dose of D = 0.1 mmol/kg body weight Gd-DTPA was injected at the start of the fifth acquisition using a power injector. This study was provided by the Paul Strickland Scanner Centre at Mount Vernon Hospital, Northwood, UK. Data from this study was acquired in accordance with the recommendation given by Leach et al. [32] and previously reported [33], [14]. Tumor ROIs were drawn by an expert radiologist based on the dynamic T 1 images.

III. RESULTS
To validate the proposed methods, the first simulation study evaluates the fit of the semi-parametric technique to simulated data in comparison to a parametric method. The second simulation study evaluates the estimation algorithm for the onset of the enhancement and the kinetic parameters when the arrival of tracer in the tissue ROI is lagged. For in-vivo validation, clinical data of 24 DCE-MRI scans of breast cancer patients were analyzed. A series of dynamic images from a patient is depicted in Fig. 1.

A. Simulation studies
The semi-parametric technique provides an accurate fit to the observed contrast agent concentration time curve.
The sum of the squared residuals (SSR) for the simulated data has a range of 1.89 · 10 −4 to 2.13 · 10 −3 with a mean of 1.13 · 10 −3 over all the 13 experiments. With the reference parametric model, the SSR has a range of 0.198 to 1.520 with a mean of 0.661, i.e., the fit to the observed data was poor in the reference parametric model compared to the proposed semi-parametric technique. Fig. 2 depicts the the kinetic parameter estimates for the semi-parametric technique and the reference parametric method compared to the ground truth. Fig. 2 (a) shows that the parameter K trans is underestimated with the reference parametric model -on average by 17.7%, which is consistent with previously published results [6]. By contrast, K trans estimates with the semi-parametric technique are much closer to the ground truth. The mean deviation from the ground truth is 6.2%. These results suggest that the semi-parametric technique is more stable compared to the parametric model.changes in v p .
With the proposed semi-parametric technique, the parameters k ep and v p are also estimated accurately, but k ep is slightly overestimated by 1.5% to 6.2% compared to a strong overestimation of 39% to 100% with the reference parametric method. For most experiments, the v p parameter obtained from the semi-parametric technique is much closer to the ground truth than that of the the reference model.  The MAD is stable with increasing time lag, but changes periodically with respect to the sampling interval of 8 seconds. For comparison, the lagged data was further analyzed without t 0 estimation. For a small time lag up to 8 seconds, i.e., up to one sampling interval, a lag of the onset time has negligible influence on K trans estimation.
However, for larger onset lags MAD increases significantly when onset is not taken into account. Table III shows the sum of squared residuals (SSR), i.e., the goodness of fit for the data from all 24 in-vivo scans obtained from the semi-parametric and the parametric approach. As with the simulated data, the fit is clearly better for the semi-parametric technique. Table IV  In this study, the 95% CI of the contrast agent concentration time series nicely covers the observed CTC. About half the observed points, however, lie outside the CI. The estimated observation error is larger than the estimated error of the CTC, i.e., the uncertainty of the model about the true CTC. For the first scan of subject 4, the CI is noticeably wider than that of other scans, i.e., the quality of the data from this scan is lower than the quality of the other scan. Typically, the CI is wider after the upslope of the CTC, and narrow during wash-out. It widens again at the end of the time series, because the spline fit at the end point relies on less data. the estimated response always has the expected shape: start at zero, rapid upslope, a short plateau and a nearly exponential downslope. The response function is well characterized by using Eqn. 20, and gives robust estimates of the kinetic parameters across the given ROI. However, the first part of the response function is typically slightly underestimated and the second part is slightly overestimated. A more flexible parametric model might fit the response function better. Fig. 5 (rows 1 and 2) depicts the sum of squared residuals (SSR) in the ROI for the semi-parametric technique and the reference parametric method for the pre-treatment scans of three patients. SSRs are clearly reduced over the whole area of the tumor, especially in areas with high K trans values, i.e., in areas with fast blood flow. Therefore, K trans estimates are clearly different for both methods, parameter maps are shown in Fig. 5 (row 3 and 4). There is, however, no general trend for under-or overestimation with different methods. For example, for subject 9, K trans estimates are higher with the parametric method, as the upslope is overestimated (see also Fig. 4). On the other hand, for patient 3, K trans estimates are to low by using the parametric technique due to a wrong onset time. Again, the Bayesian method allows us to access the precision of the estimated kinetic parameters. Non-linear parametric models are typically difficult to estimate due to the convergence issues and problems in specifying consistent starting values. Bayesian non-linear regression can overcome these convergence problems, but at a cost of computational time [14]. The proposed semi-parametric method provides a reliable and practical way of circumventing these problems. With the proposed technique, computation only includes sampling from standard distributions (see section II-B.4) which can be done efficiently [35] and gives good mixing of the MCMC kernel.

B. In vivo validation
In contrast to classical approaches [12], Bayesian P-splines allow simultaneous estimation of model and smoothing parameters. Adaptive smoothing parameters can be obtained and are important to model the sharp changes in the dynamic series of the contrast concentration.
One unique feature of the paper is the estimation of the onset of the CTC based on Bayesian inference. Onset time can be determined for the whole region of interest or on a more local level. It is an important parameter in quantitative DCE-MRI models as it has great influence on other parameters in the kinetic model and an incorrect onset time can lead to a strong bias in parameter estimates [15]. Onset time is also an important clinical index for characterizing suspicious breast lesions [36], although in this study they are not evident for the patients studied.
The proposed method allows a direct assessment of the response function, i.e., the actual flow in the tissue. can therefore serve as a supporting tool for quality control. When choosing a pixel, the CTC would be displayed interactively along with the semi-parametric fit, the CI, and the estimated onset time. This would not only help to reassess the quality of the data, but also make it easier to understand what is actually happening physiologically.