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

IEEE TRANSA CTIONS ON MEDICAL IMA GING 1 A semi-parametric technique for the quantitati v e analysis of Dynamic contrast-enhanced MR images based on Bayesian P-Splines V olker J Schmid, Brandon Whitcher , Anwar R Padhani and Guang-Zhong Y ang ∗ Abstract Dynamic Contrast-enhanced Magnetic Resonance Imaging (DCE-MRI) is an important tool for detecting subtle kinetic changes in cancerous tissue. Quantitativ e analysis of DCE-MRI typically in volv es the con volution 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, con vergence of the optimization algorithm is not guaranteed and the accuracy of the model fitting may be compromised. T o ov ercome this problems, this paper proposes a semi-parametric penalized spline smoothing approach, with which the AIF is conv olved with a set of B-splines to produce a design matrix using locally adaptiv e smoothing parameters based on Bayesian penalized spline models (P-Splines). It has been shown that kinetic parameter estimation can be obtained from the resulting decon volv ed response function, which also includes the onset of contrast enhancement. Detailed validation of the method, both with simulated and in vivo data, is provided. Index T erms Bayesian hierarchical modeling, dynamic contrast-enhanced magnetic resonance imaging, onset time, penalty splines, pharmacokinetic models, semi-parametric models I . I N T RO D U C T I O N Ev aluation of tissue kinetics in cancer with Dynamic Contrast-Enhanced Magnetic Resonance Imaging (DCE- MRI) has become an important tool for cancer diagnosis and quantification of the outcome of cancer therapies V olker Schmid and Guang-Zhong Y ang* are with the Institute of Biomedical Engineering, Imperial College, London SW7 2AZ, United Kingdom. Brandon Whitcher is with the Clinical Imaging Centre, GlaxoSmithKline, London, United Kingdom. Anwar Padhani is with Paul Strickland Scanner Centre, Mount V ernon Hospital, Northwood, United Kingdom. Asterisk indicates corresponding author (e-mail g.z.yang@imperial.ac.uk). Preliminary results initially presented at 9th International Conference on Medical Image Computing and Computer Assisted Intervention, October 1st to 6th 2006, Copenhagen [1]. 2 IEEE TRANSA CTIONS ON MEDICAL IMA GING [2]. For DCE-MRI, after the administration of a contrast agent, such as Gadolinium diethylenetriaminepentaacetic acid (Gd-DTP A), a dynamic imaging series is acquired. T ypically , 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]. W ith Gd-DTP A, the agent does not enter into cells, so DCE-MRI depicts exchange between the vascular space and the extra-v ascular extra-cellular space (EES). Current approaches to quantitati ve analysis of DCE-MRI generally rely on non-linear pharmacokinetic models. These models are usually deri ved from the solution to a system of linear differential equations, which describe the blood flo w in the tissue[5]. In practice, single-compartment models may not always be suitable, and thus more complex models ha ve been suggested. They include the “e xtended” T ofts–Kermode model (see Eqn. 2) and the tissue homogeneity approach [6], [7]. Howe ver , 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 hav e recei ved extensi ve attention in quantitati ve 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 model- free approach. The formulation of the response function based on a B-spline polynomial representation, howe ver , is ill-conditioned and a first-order difference penalty spline (Tikhono v regularization) has been imposed. As an alternati ve, Bayesian inference has also been in vestigated to replace traditional least-square fitting algo- rithms in MRI [13], [14]. Bayesian methods allo w a more accurate description of the estimation uncertainty and can ef fecti vely 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 relativ ely 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 v alidation can be used. In Bayesian frame works, P-splines regression parameters can be estimated jointly with the penalty weight and hence allow for adapti ve 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 adapti ve smoothing approach as the observ ed time series signal v aries rapidly in the first minute after injecting the contrast bolus. The proposed algorithm provides a semi-parametric decon volution approach and results in a smooth response function, along with the corresponding estimates of uncertainty . Kinetic parameters are rob ustly deriv ed 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 deriv ed from the Bayesian estimation process. Detailed v alidation of the method both with simulated and in vivo data of patients with breast tumors is provided. A SEMI-P ARAMETRIC TECHNIQUE FOR THE Q U ANTIT A TIVE AN AL YSIS OF DCE-MRI B ASED ON BA YESIAN P-SPLINES 3 I I . T H E O RY A N D M E T H O D S 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 con volution of the arterial input function (AIF) with a single exponential function [23], i.e. , C t ( t ) = C p ( t ) ⊗ K trans exp( − k ep t ) . (1) 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 T ofts–K ermode model [6] as a reference parametric model gi ven by C t ( t ) = v p C p ( t ) + C p ( t ) ⊗ K trans exp( − k ep t ) , (2) 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 giv en by [24] C p ( t ) = D 2 X i =1 a i exp( − m i t ) (3) 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. W ith this explicit form of the AIF , the con volution in Eqn. 2 can be deriv ed analytically , i.e. , C t ( t ) = v p D 2 X i =1 a i exp( − m i t ) + D K trans 2 X i =1 a i { exp( − m i t ) − exp[ − k ep t ] } k ep − m i . (4) The kinetic parameters K trans , k ep and v p are estimated by fitting Eqn. 4 to the observed data. Optimization is usually performed using the Lev enberg–Marquardt algorithm [26]. B. Bayesian P-Spline model for DCE-MRI Parametric models, howe ver , may not giv e an accurate fit of the observed data, as they are too simplistic and can ov erlook ef fects such as flo w heterogeneity or the water exchange ef fect [27], [28]. More complex models have since been de veloped [7], but they are more difficult to optimize. It has been found that semi-parametric models can fit the data more accurately and in this study penalty splines (P-Splines) in a Bayesian hierarchical framework are used. 4 IEEE TRANSA CTIONS ON MEDICAL IMA GING 1) Discrete Decon volution: Mathematically , a more general expression of Eqn. 1 can be written as: C t ( t ) = C p ( t ) ⊗ f ( t ) , (5) where f ( t ) is the response function in the tissue. Assuming that C p ( · ) and f ( · ) are constant o ver small interv als ∆ t , a discretized form of Eqn. 5 is giv en by C t ( τ i ) = T X j =1 C p ( τ i − t j ) f ( t j )∆ t = T X j =1 A ij f ( t j ) , (6) where C t is measured on discrete time points τ 1 . . . , τ n . The n × T matrix A may be interpreted as a con volution operator and is defined by A ij =    C p ( t n i − j +1 )∆ t if τ i ≤ t j ; 0 otherwise , (7) 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 e v aluated from Eqn. 3 – at time points t 1 , . . . , t T , which can be different from τ 1 , . . . , τ n . 2) P enalty Splines: By solving Eqn. 6, the response function f ( t ) can be deconv olved from C p ( t ) . Ho wev er, this system may be numerically unstable, i.e. , the decon volved response function is susceptible to noise. T o ov ercome this problem, we assume that f ( t ) is a smooth, k -times dif ferentiable function. T o this end, a B-spline representation of the response function is used in this study , i.e. , f ( t ) = p X j =1 β j B tj , (8) where B is the n × p design matrix of k th order B-splines with knots s 1 , . . . , s p + k . In vector notation f = ( f ( t 1 ) , . . . , f ( t T )) 0 and Eqn. 8 may be expressed as f = B β . (9) Accordingly , Eqn. 6 can be written as C t = Af = AB β = D β , (10) where D = AB is a n × p design matrix, representing the (discrete) conv olution of the AIF with the B-spline polynomials. T o enhance the numerical stability , a penalty on the B-spline regression parameters is introduced, such that β t = 2 β t − 1 − β t − 2 + e t for t = 3 , . . . , p. (11) 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 adaptiv e. T o this end, we use a Bayesian hierarchical framew ork for parameter inference. A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 5 3) Bayesian hierar chical framework: In Bayesian inference, a priori information, i.e. , information av ailable before observation of measurement, has to be expressed in terms of probability distributions. Here, our prior kno wledge is the assumption about the model. In this paper , we assume that the observ ed contrast agent concentration C t is noisy realization of the true model (Eqn. 10), i.e. , C t ( t ) ∼ N ( D t β , σ 2 ) for all t, (12) where D t denotes the t th ro w of D . That is, a priori the error is assumed to be Gaussian distributed with unknown v ariance σ 2 . W e use a relatively flat prior for the variance parameter , i.e. , σ 2 ∼ IG (1 , 10 − 5 ) , (13) where IG denotes the In verse Gamma distrib ution. The penalty in Eqn. 11 can be expressed as a priori distribution on β [21], β t ∼ N (2 β t − 1 − β t − 2 , δ 2 t ) for t = 3 , . . . , p, (14) where δ 2 t is the variance of e t . For a locally adaptive estimation of δ 2 t , the following prior model is used δ 2 i ∼ IG (10 − 5 , 10 − 5 ) for i = 3 , . . . , p, (15) which allows time v arying 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, giv en by p ( β , δ 2 , σ 2 ) ∝ ` ( C t | β , σ 2 ) p ( σ 2 ) p ( β | δ 2 ) p ( δ 2 ) . (16) 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 distrib ution [29]. The full conditional of β , i.e. , the joint distribution of β gi ven all other parameters and the data, C t , is a p -dimensional multiv ariate normal distribution β | C t , δ 2 , σ 2 ∼ N p  C t 0 D ( D 0 D + R ) − 1 , σ 2 ( D 0 D + R ) − 1  , (17) where R is the in verse cov ariance matrix of the prior distribution of β [21]. The full conditionals of both variance parameters are independent In verse Gamma distributions δ 2 t | β iid ∼ IG  10 − 5 + 0 . 5 , 10 − 5 + 0 . 5( β t − 2 β t − 1 + β t − 2 )  for t = 3 , . . . , p, (18) σ 2 | β , C t ∼ IG  1 + T 2 , 10 − 5 + 0 . 5( C t − D t β ) 0 ( C t − D t β )  . (19) T o assess the posterior distribution, samples of the parameters β , σ 2 and δ 2 are drawn alternately from Eqn. 17, Eqn. 18 and Eqn. 19. After a suf ficient burn-in period, the MCMC algorithm produces samples from the posterior distribution. 6 IEEE TRANSA CTIONS ON MEDICAL IMA GING C. Estimating kinetic parameters fr om 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-con volved 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. W e make use of the fact that the Bayesian approach pro vides 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 arri val of the tracer in the tissue of interest is unknown. Howe ver , 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 cov er 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 deriv ativ e of C t ( t ∗ ) can be computed by d dt C t ( t ) = d dt [ C p ( t ) ⊗ f ( t )] = C p ( t ) ⊗ d dt f ( t ) . Since f ( t ) is a spline, the deri v ativ e can be computed as [31] f ( t ) = n − k X j =1 γ j B ( k − 1) t ( j +1) , where B ( k − 1) is the design matrix of ( k − 1) th order B-Splines and γ j is defined via γ j = k s j + k +1 − s j +1 ( β j +1 − β j ) for j = 1 , . . . , J − 1 . 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 t 0 = t ∗ − C t ( t ∗ ) dC t ( t ∗ ) /dt . DCE-MRI studies typically assume that the onset of the enhancement is the same ov er the whole region of interest (R OI). In this case, the median of the estimated t 0 for all voxels in the R OI may be used as an estimate of the global v alue. Howe ver , for larger R OIs, local estimates of the onset may be required, and can be computed from the proposed technique. A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 7 2) Obtaining kinetic parameters: The semi-parametric technique provides a de-con volv ed and de-noised response function. In order to obtain kinetic parameters, we can fit a non-linear function to the estimated response. T o this end, the following model has been used: f ( t ) = F p ·          E exp[ − ( t − t 0 − T c ) E F p /v e ] for t ≥ ( T c + t 0 ) , 1 for t 0 ≤ t < ( T c + t 0 ) , 0 for t < t 0 , (20) where T c = v p /K trans is the transit time through the capillary , v e = K trans /k ep is the v olume 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 (AA TH) [7]. In this model extraction fraction E and mean plasma flow F p may be not identifiable, but the product E F 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 Lev enberg–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]. W e use the first series which was designed to be representativ e of data acquired from a breast tumor . Data was simulated using MMID4, part of a softw are made a v ailable by the National Simulation Resource, Dept. of Bioengineering, Uni versity of W ashington ( http://www .nsr .bioeng .washington.edu ). Estimates from the literature were used as baseline values. For twelve further experiments, one of the kinetic parameters F p , v p and P S was changed four times while the other two where held fixed at the baseline values (see T ab . I). Data was sampled at 1 Hz. For a second set of simulated data, a time lag of up to 30 seconds was added to the simulated data to e v aluate the ef fect of lagged contrast uptake. T o make the experiment more realistic for DCE-MRI, data was down sampled to 1/8 Hz; scans in DCE-MRI experiments are typically acquired ev ery 4–12 seconds. In vivo data was deriv ed 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 MA GNETOM 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-DTP A 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 V ernon Hospital, Northwood, UK. Data from this study was acquired in accordance with the recommendation gi ven by Leach et al.[32] and previously reported [33], [14]. T umor R OIs were drawn by an expert radiologist based on the dynamic T 1 images. 8 IEEE TRANSA CTIONS ON MEDICAL IMA GING I I I . R E S U L T S T o v alidate the proposed methods, the first simulation study ev aluates the fit of the semi-parametric technique to simulated data in comparison to a parametric method. The second simulation study ev aluates the estimation algorithm for the onset of the enhancement and the kinetic parameters when the arriv al 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 ov er all the 13 experiments. W ith 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) sho ws 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 . W ith the proposed semi-parametric technique, the parameters k ep and v p are also estimated accurately , but k ep is slightly ov erestimated by 1.5% to 6.2% compared to a strong ov erestimation 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. For small values of v p (experiment 6), howe ver , the semi-parametric method shows a larger deviation from ground truth. An important advantage of the proposed Bayesian technique is that it not only produces point estimates, b ut also estimations of the errors or interv al estimators. T ab. II gi ves 95% Credible Intervals (CI) for K trans for all the 13 experiments conducted. For experiments 8 and 9, where v p has larger values, the 95% CI are broad, but still co ver the true value. Similar CIs are av ailable for the other parameters. The Bayesian technique provides important information about both the accuracy and precision of its estimates. T o v alidate the proposed algorithm for estimating the onset of the contrast uptake, a time lag of up to 30 seconds was added to the down-sampled simulated data. Data was analyzed with the proposed semi-parametric technique and the enhancement onset time and the kinetic parameters were deriv ed by equations described in section II-C. The estimation of t 0 with the proposed method is sho wn to be accurate. Correlation between the estimated onset time and ground truth is 0 . 9984 . The mean difference between the true and estimated onset time is 0 . 1800 seconds with a standard de viation of 0 . 9339 seconds and a maximum absolute difference of 2 . 2992 seconds. A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 9 Fig. 3 shows the Mean Absolute Dif ference (MAD) between the estimated v alues of K trans and the ground truth. The MAD is stable with increasing time lag, but changes periodically with respect to the sampling interv al 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. Ho wev er, for larger onset lags MAD increases significantly when onset is not taken into account. B. In vivo validation T able 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. T able IV provides the estimated enhancement onset time for all scans; the scanning protocol indicates that the contrast agent was injected at the start of the fifth image, i.e. after 49 . 40 s . For three subjects, the onset time is displayed in Fig. 4 as vertical dashed line. V isual inspection of this figure and of the other subjects shows that the estimated onset time is consistent with the observed contrast concentration time series. F or most scans, the estimated onset time is between the start of the fifth acquisition and the start of the sixth acquisition, which is reasonable giv en the scanning protocol. For some scans, ho wever , the onset time is much larger due to deviation from the scanning protocol. Figs. 4 depicts the observed contrast concentration time series for three voxels from three different scans together with the estimates from the parametric reference method and those from the proposed semi-parametric approach with 95% CI. It is e vident that in all vox els the semi-parametric method fits the data better than the parametric approach. The figures also indicate that the standard parametric model is not always appropriate for the observed data, where the initial upslope is difficult to follow by using the parametric model. The model also fails when the assumed onset time strongly deviates from the observed onset time, for example for the first scan of subject 4 (see bottom row of Fig. 4). In this study , the 95% CI of the contrast agent concentration time series nicely covers the observed CTC. About half the observed points, howe ver , 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. T ypically , 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. Fig. 4 also depicts the median of the estimated de-con volved response function and the fit to the parametric model Eqn. 20 for the same three v oxels. Although there is no restriction on the form of the response function, 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 gi ven ROI. Howe ver , the first part of the response function is typically slightly underestimated and the second part is slightly ov erestimated. A more fle xible parametric model might fit the response 10 IEEE TRANSA CTIONS ON MEDICAL IMA GING function better . Fig. 5 (rows 1 and 2) depicts the sum of squared residuals (SSR) in the R OI for the semi-parametric technique and the reference parametric method for the pre-treatment scans of three patients. SSRs are clearly reduced ov er the whole area of the tumor , especially in areas with high K trans v alues, 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, ho wev er, 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 allo ws us to access the precision of the estimated kinetic parameters. Fig. 5 (ro w 5) depicts the standard error of the K trans parameter computed with the semi-parametric technique. Estimation errors are higher with higher values, but they are generally relativ ely small. For most vox els, the relative error, i.e. , standard error di vided by estimated parameter , does not exceed 20%. I V . C O N C L U S I O N In this paper, we hav e introduced a semi-parametric technique based on Bayesian P-spline models for quantifying CTCs obtained with DCE-MRI. Compared with parametric models, the proposed semi-parametric technique pro vides a superior fit to the observed concentration time curves. In particular, the proposed semi-parametric method captures the upslope of the time series accurately , which is important for an accurate fit of the CTCs in DCE-MRI and proper calculation of K trans . Results from the simulated data sho w a clear improvement in both time curve fit and parameter estimation. In addition, the Bayesian method provides information about the precision of the estimation, including standard errors or credible intervals. For in vivo v alidation, the fit of the contrast agent concentration with the proposed semi- parametric technique is superior to the parametric method. Estimates of kinetic parameters for both techniques are dif ferent where the fit of the parametric technique is poor , suggesting that the estimation of kinetic parameters with the proposed semi-parametric technique is more accurate in these areas. The exact assessment of the kinetics is clinically important as changes in respect to the response function and the kinetic parameters can be subtle, especially for drugs that cause multiple effect on tumor vasculature, such as combinations of antiangiogenic drugs [34]. Non-linear parametric models are typically difficult to estimate due to the con ver gence issues and problems in specifying consistent starting v alues. Bayesian non-linear regression can overcome these con vergence problems, b ut at a cost of computational time [14]. The proposed semi-parametric method provides a reliable and practical way of circumventing these problems. W ith 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. In contrast to classical approaches [12], Bayesian P-splines allow simultaneous estimation of model and smoothing parameters. Adapti ve smoothing parameters can be obtained and are important to model the sharp changes in the A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 11 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 quantitati ve 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. Parameters of interest may be estimated by fitting a non-linear model to the response function. Smoothing via P-splines provides an effecti ve way of error reduction and decon volution of the arterial input function. The proposed technique also allows the quantification of errors both in fitting the observed data and in estimating kinetic parameters. Thus far , the estimation error in DCE-MRI models is rarely discussed, the proposed technique can contribute to the e valuation of the quality of DCE-MRI scans. Results from the semi-parametric technique can therefore serve as a supporting tool for quality control. When choosing a pixel, the CTC would be displayed interacti vely 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, b ut also make it easier to understand what is actually happening physiologically . A C K N O W L E D G M E N T S Support for V olker Schmid was supported by a platform grant (GR/T06735/01) from the UK Engineering and Physical Sciences Research Council (EPSRC). He was also partially financed through a research grant from GlaxoSmithKline. W e are grateful to David Buckley at Imaging Science and Biomedical Engineering, Uni versity of Manchester , UK for providing the simulated data. The clinical data were graciously provided by Anwar Padhani and Jane T aylor at Paul Strickland Scanner Centre, Mount V ernon Hospital, Northwood, UK. R E F E R E N C E S [1] V . J. Schmid, B. Whitcher, and G. Z. Y ang, “Semi-parametric analysis of dynamic contrast-enhanced MRI using Bayesian P-splines, ” in Medical Imaging Computing and Computer-Assisted Intervention – MICCAI 2006 , ser . Lecture Notes in Computer Science, R. Larsen, M. Nielsen, and J. Sporring, Eds., no. 4190. Berlin: Springer , 2006, pp. 679–686. [2] D. J. Collins and A. R. Padhani, “Dynamic magnetic resonance imaging of tumor perfusion, ” IEEE Engineering in Biology and Medicine , vol. 23, no. 5, pp. 65–83, Sep-Oct 2004. [3] G. J. M. Parker and A. R. Padhani, “ T 1 -w DCE-MRI: T 1 -weighted dynamic contrast-enhanced MRI, ” in Quantitative MRI of the Brain , P . T ofts, Ed. Chichester , England: Wile y , 2003, ch. 10, pp. 341–364. [4] E. K. Fram, R. J. Herfkens, G. A. Johnson, G. H. Glover , J. P . Karis, A. Shimakawa, T . G. Perkins, and N. J. Pelc, “Rapid calculation of T1 using variable flip angle gradient refocused imaging, ” Magnetic Resonance Imaging , v ol. 5, pp. 201–208, 1987. [5] H. B. Larsson and P . S. T ofts, “Measurement of the blood-brain barrier permeability and leakage space using dynamic Gd-DTP A scanning – a comparison of methods, ” Magnetic Resonance in Medicine , vol. 24, no. 1, pp. 174–176, 1992. [6] D. L. Buckley , “Uncertainty in the analysis of tracer kinetics using dynamic contrast-enhanced T 1 -weighted MRI, ” Magnetic Resonance in Medicine , vol. 47, no. 3, pp. 601–606, 2002. 12 IEEE TRANSA CTIONS ON MEDICAL IMA GING [7] K. S. St.Lawrence and T .-Y . Lee, “ An adiabatic approximation to the tissue homogeneity model for water exchange in the brain: I. theoretical deriv ation, ” Journal of Cerebral Blood Flow & Metabolism , vol. 18, pp. 1365–1377, 1998. [8] T . T wellmann, O. Lichte, and T . W . Nattkemper, “ An adaptiv e tissue characterization network for model-free visualization of dynamic contrast-enhanced magnetic resonance image data, ” IEEE T ransactions on Medical Imaging , vol. 24, no. 10, pp. 1256–66, October 2005. [9] K.-H. Chaung, M.-T . W u, Y .-R. Lin, K.-S. Hsieh, M.-L. W u, S.-Y . Tsai, C.-W . K o, and H.-W . Chung, “ Application of model-free analysis in the MR assessment of pulmonary perfusion dynamics, ” Magnetic Resonance in Medicine , vol. 54, no. 2, pp. 299–308, 2005. [10] V . J. Cunningham and T . Jones, “Spectral analysis of dynamic PET studies, ” Journal of Cer ebral Blood Flow and Metabolism , vol. 13, pp. 15–23, 1993. [11] R. N. Gunn, S. R. Gunn, F . E. Turkheimer , J. A. D. Aston, and V . J. Cunningham, “Positron emission tomography compartmental models: A basis pursuit strategy for kinetic modeling, ” Journal of Cer ebral Blood Flow and Metabolism , vol. 22, pp. 1425–1439, 2002. [12] M. Jerosch-Herold, C. Swingen, and R. Seethamraju, “Myocardial blood flow quantification with MRI by model-independent decon volution, ” Medical Physics , vol. 29, no. 5, pp. 886–897, May 2002. [13] C. G ¨ ossl, L. Fahrmeir , and D. P . Auer , “Bayesian modeling of the hemodynamic response function in BOLD fMRI, ” NeuroIma ge , vol. 14, pp. 140–148, 2001. [14] V . J. Schmid, B. Whitcher , G. Z. Y ang, N. J. T aylor, and A. R. Padhani, “Statistical analysis of pharmacokinetical models in dynamic contrast-enhanced magnetic resonance imaging, ” in Medical Imaging Computing and Computer -Assisted Intervention – MICCAI 2005 , ser . Lecture Notes in Computer Science, J. Duncan and G. Gerig, Eds., no. 3750. Berlin: Springer , 2005, pp. 886–893. [15] M. R. Orton, D. J. Collins, S. W alker -Samuel, J. A. d’Arcy , D. J. Hawkes, D. Atkinson, and M. O Leach, “Bayesian estimation of pharmacokinetic parameters for DCE-MRI with a robust treatment of enhancement onset time, ” Physics in Medicine and Biology , vol. 52, pp. 2393–2408, 2007. [16] V . J. Schmid, B. Whitcher, A. R. Padhani, N. J. T aylor , and G.-Z. Y ang, “Bayesian methods for pharmacokinetic models in dynamic contrast-enhanced magnetic resonance imaging, ” IEEE T ransactions on Medical Imaging , vol. 25, no. 12, pp. 1627–1636, December 2006. [17] P . Eilers and B. Marx, “Flexible smoothing with B-splines and penalties (with comments and rejoinder), ” Statistical Science , vol. 11, no. 2, pp. 89–121, 1996. [18] B. Marx and P . Eilers, “Direct generalized additi ve modeling with penalized likelihood, ” Computational Statistics and Data Analysis , vol. 28, no. 2, pp. 193–209, 1998. [19] P . R. Johnston and R. M. Gulrajani, “Selecting the corner in the L-curve approach to Tikhonov regularization, ” IEEE T ransaction on Biomedical Engineering , vol. 47, p. 1293, Sept. 2000. [20] S. Lang, E. Fronk, and L. Fahrmeir , “Function estimation with locally adaptiv e dynamic models, ” Computational Statistics , vol. 17, pp. 479–500, 2002. [21] S. Lang and A. Brezger, “Bayesian P-splines, ” Journal of Computational and Graphical Statistics , vol. 13, pp. 183–212, 2004. [22] S. Kety , “Blood-tissue exchange methods. Theory of blood-tissue exchange and its applications to measurement of blood flow . ” Methods in Medical Resear ch , vol. 8, pp. 223–227, 1960. [23] P . S. T ofts, G. Brix, D. L. Buckley , J. L. Evelhoch, E. Henderson, M. V . Knopp, H. B. W . Larsson, T .-Y . Lee, N. A. Mayr , G. J. M. Parker , R. E. Port, J. T aylor , and R. W eiskoff, “Estimating kinetic parameters from dynamic contrast-enhanced T 1 -weighted MRI of a diffusable tracer: Standardized quantities and symbols, ” Journal of Magnetic Resonance Imaging , vol. 10, pp. 223–232, 1999. [24] P . T ofts and A. Kermode, “Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging – 1. Fundamental concepts, ” Magnetic Resonance in Medicine , vol. 17, pp. 357–367, 1991. [25] T . Fritz-Hansen, E. Rostrup, K. B. Larsson, L. Søndergaard, P . Ring, and O. Henriksen, “Measurement of the arterial concentration of Gd-DTPA using MRI: A step tow ard quantitativ e perfusion imaging, ” Magnetic Resonance in Medicine , vol. 36, pp. 225–231, 1996. [26] J. J. Mor ´ e, “The Le venber g-Marquardt algorithm: Implementation and theory , ” in Numerical Analysis: Pr oceedings of the Biennial A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 13 Confer ence held at Dundee, J une 28-J uly 1, 1977 (Lectur e Notes in Mathematics #630) , G. A. W atson, Ed. Berlin: Springer-V erlag, 1978, pp. 104–116. [27] K. Kroll, N. W ilke, M. Jerosch-Herold, Y . W ang, Y . Zhang, R. J. Bache, and J. B. Gassingthwaighte, “Modeling regional myocardial flows from residue functions of an intravascular indicator , ” American J ournal of Heart Circulation Physiology , vol. 271, pp. 1643–1655, 1996. [28] C. Landis, X. Li, F . W . T elang, J. A. Coderre, P . L. Micca, W . D. Rooney , L. L. Latour , G. V ´ etek, I. P ´ alyka, and C. S. Springer , “Determination of the MRI contrast agent concentration time course in vi vo following bolus injection: Effect of equilibrium transcytolemmal water exchange, ” Magnetic Resonance in Medicine , vol. 44, pp. 567–574, 2000. [29] W . R. Gilks, S. Richardson, and D. J. Spiegelhalter , Eds., Markov Chain Monte Carlo in Practice . Chapman & Hall, London, 1996. [30] M. Jerosch-Herold, X. Hu, N. Murthy , and R. Seethamraju, “T ime delay for arriv al of MR contrast agent in collateral-dependent myocardium, ” IEEE T ransactions on Medical Imaging , vol. 23, no. 7, pp. 881–890, July 2004. [31] N.-K. Tsao and T .-C. Sun, “On the numerical computation of the derivati ves of a B-spline series, ” IMA Journal of Numerical Analysis , vol. 13, no. 3, pp. 343–364, 1993. [32] M. O. Leach, K. M. Brindle, J. L. Evelhoch, J. R. Griffiths, M. R. Horsman, A. Jackson, G. C. Jayson, I. R. Judson, M. V . Knopp, R. J. Maxwell, D. McIntyre, A. R. Padhani, P . Price, R. Rathbone, G. J. Rustin, P . S. T ofts, G. M. T ozer, W . V ennart, J. C. W aterton, S. R. W illiams, and P . W orkman, “The assessment of antiangiogenic and antiv ascular therapies in early-stage clinical trials using magnetic resonance imaging: issues and recommendations, ” British J ournal of Cancer , vol. 92, pp. 1599–1610, 2005. [33] M.-L. W . Ah-See, A. Makris, N. J. T aylor, R. J. Burcombe, M. Harrison, J. J. Stirling, P . I. Richman, M. O. Leach, and A. R. Padhani, “Does vascular imaging with MRI predict response to neoadjuvant chemotherapy in primary breast cancer?” J ournal of Clinical Oncology (Meeting Abstracts) , vol. 22, no. 14S, p. 582, 2004. [34] A. R. Padhani and M. O. Leach, “ Anti vascular cancer treatments: functional assessments by dynamic contrast-enhanced magnetic resonance imaging, ” Abdominal Imaging , vol. 30, no. 3, pp. 325–342, 2005. [35] H. Rue, “Fast sampling of Gaussian Markov random fields, ” Journal of the Royal Statistical Society Series B , vol. 63, pp. 325–338, 2001. [36] C. Boetes, J. Barentsz, R. Mus, R. v an der Sluis, L. van Erning, J. Hendriks, R. Holland, and S. Ruys, “MR characterization of suspicious breast lesions with a gadolinium-enhanced T urboFLASH subtraction technique, ” Radiology , vol. 193, no. 3, pp. 777–781, 1994. 14 IEEE TRANSA CTIONS ON MEDICAL IMA GING T ABLE I V A L U E S O F T H E K I N E T I C P A R A M E T E R S D E S C R I B I N G T H E B E H A V I O R O F M M I D 4 U S E D T O S I M U L AT E 1 3 D I FF E R E N T C O N T R A S T C O N C E N T R ATI O N T I M E C U RV E S R E P R E S E N TA T I V E O F A B R E A S T T U M O R . Exp. Baseline 2 3 4 5 6 7 8 9 10 11 12 13 F p 0.57 0.17 0.37 0.77 0.97 0.57 0.57 v p 0.06 0.06 10 − 4 0.03 0.09 0.12 0.06 P S 0.33 0.33 0.33 0.01 0.17 0.49 0.65 v e 0.45 0.45 0.45 0.45 T ABLE II T R U E A N D E S T I M AT E D V A L U E S O F K T R A N S W I T H 9 5 % C I . Experiment 1 2 3 4 5 6 7 8 9 10 11 12 13 T rue K trans 0.251 0.146 0.218 0.268 0.280 0.233 0.233 0.233 0.233 0.010 0.147 0.323 0.388 2.5% quantile K trans 0.208 0.128 0.201 0.236 0.242 0.181 0.220 0.224 0.229 0.004 0.129 0.293 0.347 Median posterior K trans 0.245 0.138 0.218 0.263 0.284 0.200 0.237 0.281 0.282 0.005 0.147 0.320 0.376 97.5% quantile K trans 0.267 0.167 0.240 0.309 0.376 0.239 0.256 0.334 0.372 0.060 0.183 0.356 0.415 T ABLE III S U M O F S Q U A R E D R E S I D U A L S ( S S R ) B E T W E E N DA TA A N D FI T T E D F U N C T I O N F O R PA R A M E T R I C A N D S E M I - PA R A M E T R I C M E T H O D , A V E R AG E D OV E R A L L VO X E L S I N T H E RO I . 1st scan, subject 1 2 3 4 5 6 parametric 0.0682 0.1579 0.2209 0.1825 0.6101 0.4382 semi-parametric 0.0217 0.0575 0.0341 0.183 0.0398 0.0333 7 8 9 10 11 12 parametric 0.0723 0.0366 0.0798 0.0568 0.1420 0.0552 semi-parametric 0.0242 0.0068 0.104 0.0077 0.0131 0.0234 2nd scan, subject 1 2 3 4 5 6 parametric 0.0457 0.0691 0.0403 0.8665 0.6469 0.5804 semi-parametric 0.0123 0.0328 0.0122 0.0484 0.0386 0.0334 7 8 9 10 11 12 parametric 0.0998 0.0525 0.0237 0.0544 0.0225 0.0289 semi-parametric 0.0209 0.0102 0.0071 0.0205 0.0072 0.0168 T ABLE IV E S T I M AT E D C O N T R A S T C O N C E N T R ATI O N E N H A N C E M E N T O N S E T T I M E I N S E C O N D S F O R A L L 2 4 S C A N S . P R OT O C O L S P E C I FI E D T H AT T H E T R A C E R W A S I N J E C T E D A F T E R 4 9 . 4 0 S E C O N D S . Subject 1 2 3 4 5 6 7 8 9 10 11 12 1st scan 55.19 54.14 69.09 90.20 60.04 57.63 52.56 54.95 55.85 56.72 47.48 50.34 2nd scan 55.09 56.58 87.09 56.58 56.71 55.31 53.10 55.46 53.33 53.96 54.06 49.59 A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 15 Fig. 1. T ime series of contrast concentration images at baseline and after 23, 70, 222 and 410 seconds - central slice from first scan of patient 9. true Ktrans estimated 0.0 0.1 0.2 0.3 0.4 0.5 0.0 0.1 0.2 0.3 0.4 0.5 (a) true kep estimated 0.0 0.4 0.8 0.0 0.4 0.8 (b) true vp estimated 0.00 0.04 0.08 0.12 0.00 0.04 0.08 0.12 (c) Fig. 2. Scatter plots of parameter estimates vs. true values – These figures are similar to Fig. 3 in [6]. Estimates of the parameters (a) K trans , (b) k ep and (c) v p . Bullets represent the results using the semi-parametric method, diamonds the results using the parametric technique. 0 5 10 15 20 25 30 0.00 0.05 0.10 0.15 0.20 true onset in seconds mean deviation from true Ktrans with onset estimation without onset estimation Fig. 3. Mean absolute difference between true and estimated K trans values plotted against true enhancement onset time in seconds. Solid dots are results from semi-parametric analysis without considering onset, circles are results from the proposed algorithm with onset estimation. 16 IEEE TRANSA CTIONS ON MEDICAL IMA GING 02468 0.0 0.2 0.4 0.6 0.8 time in minutes contrast concentration in mMol observed parametric semi−parametric credible interval 02468 0.0 0.1 0.2 0.3 0.4 0.5 0.6 time in minutes arbitrary units de−convolved respons e fit to response 02468 0.0 0.2 0.4 0.6 0.8 time in minutes contrast concentration in mMol observed parametric semi−parametric credible interval 02468 0.0 0.1 0.2 0.3 0.4 0.5 0.6 time in minutes arbitrary units de−convolved respons e fit to response 02468 0.0 0.2 0.4 0.6 0.8 time in minutes contrast concentration in mMol observed parametric semi−parametric credible interval 02468 0.0 0.1 0.2 0.3 0.4 0.5 0.6 time in minutes arbitrary units de−convolved respons e fit to response Fig. 4. Left column: Observed contrast agent concentration time series (dots) of a voxel in the R OI of the pre-treatment scans for three different subjects (from top to bottom: subjects 5, 9 and 4) along with fit by parametric (slashed line) and semi-parametric technique; for the semi-parametric technique the median (solid line) and the 95% CI (dotted lines) are depicted. The estimated enhancement onset time is marked as vertical line. Right column: Median estimated de-con volved response function (dots) with fit to the model Eqn. 20 (grey line) for the same vox els. A SEMI-P ARAMETRIC TECHNIQUE FOR THE QUANTIT A TIVE ANAL YSIS OF DCE-MRI BASED ON BA YESIAN P-SPLINES 17 Patient 2 Patient 3 Patient 9 SSR parametric SSR semi-parametric K trans parametric K trans semi-parametric standard error K trans Fig. 5. Parameter maps for the first scan of patients 2, 3 and 9. From top to bottom: Sum of squared residuals (SSR) map for reference parametric technique; SSR map for proposed semi-parametric technique; K trans parameter map estimated with reference parametric technique; K trans parameter map estimated with proposed semi-parametric technique standard error of K trans estimated from proposed semi-parametric technique.