Non-parametric deconvolution using Bézier curves for quantification of cerebral perfusion in dynamic susceptibility contrast MRI.

OBJECTIVE: Deconvolution is an ill-posed inverse problem that tends to yield non-physiological residue functions R(t) in dynamic susceptibility contrast magnetic resonance imaging (DSC-MRI). In this study, the use of Bézier curves is proposed for obtaining physiologically reasonable residue functions in perfusion MRI.
MATERIALS AND METHODS: Cubic Bézier curves were employed, ensuring R(0) = 1, bounded-input, bounded-output stability and a non-negative monotonically decreasing solution, resulting in 5 parameters to be optimized. Bézier deconvolution (BzD), implemented in a Bayesian framework, was tested by simulation under realistic conditions, including effects of arterial delay and dispersion. BzD was also applied to DSC-MRI data from a healthy volunteer.
RESULTS: Bézier deconvolution showed robustness to different underlying residue function shapes. Accurate perfusion estimates were observed, except for boxcar residue functions at low signal-to-noise ratio. BzD involving corrections for delay, dispersion, and delay with dispersion generally returned accurate results, except for some degree of cerebral blood flow (CBF) overestimation at low levels of each effect. Maps of mean transit time and delay were markedly different between BzD and block-circulant singular value decomposition (oSVD) deconvolution. DISCUSSION: A novel DSC-MRI deconvolution method based on Bézier curves was implemented and evaluated. BzD produced physiologically plausible impulse response, without spurious oscillations, with generally less CBF underestimation than oSVD.

Introduction

The commonly used perfusion-related parameters cerebral blood flow (CBF), cerebral blood volume (CBV) and mean transit time (MTT) are of value in the assessment of, for example, ischemic threshold level and ischemic penumbra in stroke [1, 2], as well as for characterization and grading of brain tumors [3]. In this context, dynamic susceptibility contrast magnetic resonance imaging (DSC-MRI) is a well-established technique for assessment of parameters related to cerebral perfusion and microvasculature [4, 5].

DSC-MRI is based on the injection of contrast agent (CA), with subsequent measurement of the MRI signal loss during the CA bolus passage through the tissue of interest. Quantification of perfusion parameters from DSC-MRI data is based on general tracer kinetic theory, assuming that the measured tissue concentration time curve (CTC), calculated from the signal curve, can be expressed as the convolution of the arterial input function (AIF) and the tissue residue function scaled with CBF (i.e., the tissue impulse response function) [4-6]. Inference of CBF and MTT from DSC-MRI, thus, requires accurate voxel-wise deconvolution of the observed CTC with the measured AIF. Furthermore, the shape of the tissue residue function reflects microvascular tracer retention and capillary velocity distribution [7, 8], and the residue function may contain information about flow heterogeneity and tissue oxygen extraction [9].

Deconvolution is, in general, an ill-posed inverse problem where slight alterations in the CTC may translate to a dramatic influence on the residue function [8]. Furthermore, effects on arterial bolus delay and dispersion occurring between the site of the AIF measurement and the arterial inlet of the tissue of interest need to be considered in the deconvolution procedure [10-12]. Ideally, the applied AIF should reflect the arterial input to each tissue voxel being analyzed, but the AIF is typically sampled in a major artery (such as the middle cerebral artery or its branches) to minimize partial volume effects [13]. Numerous deconvolution techniques have previously been proposed in the literature [e.g. 5, 7, 8, 14, 15]. Fourier techniques are insensitive to bolus delay but have been shown to poorly estimate high flow rates coupled with noise [5]. The algebraic, model-independent singular value decomposition (SVD) technique was shown to be a robust alternative [5], and Wu et al. later introduced the delay-insensitive block-circulant SVD (oSVD) version [15]. Apart from CBF underestimation in many practical situations, a major weakness of SVD-based deconvolution methods is the induction of non-physiological oscillations in the obtained residue function [5]. Various regularization techniques have been adopted to reduce these oscillations, for example, SVD with a fixed threshold for the singular values, oSVD incorporating an oscillation index [15] and Tikhonov regularization [16].

Mouridsen et al. proposed a vascular model (VM) which produces smooth residue functions and is believed to provide more accurate perfusion estimates than SVD-based methods [8]. The VM was implemented in a Bayesian framework modeling tracer transit times with a gamma probability distribution function. Two parameters (one for shape and another for scale) describe the distribution. Consequently, the resulting residue functions belong to a restricted family of shapes, which may impose restrictions on the applicability and versatility of the VM approach.

Mehndiratta et al. introduced a non-parametric control point interpolation (CPI) deconvolution method, implemented in a Bayesian framework [7]. With the aid of control points, CPI parallels the flexibility in residue function shape that is obtained with non-parametric techniques. Imposing constraints on the control points allows discrimination against non-physiological residue functions. However, CPI exhibits some degree of oscillatory behavior and CBF overestimation in the case of a boxcar residue function [7]. A general cause for concern is that the requirement of the residue function to pass through all control points leads to a large number of optimisation parameters (23 when using twelve control points with correction for arterial delay).

In this work, we propose a non-parametric deconvolution technique using Bézier curves, referred to as Bézier curve deconvolution (BzD). Bézier curves were independently developed by de Casteljau and Bézier in 1959–60 for generation of smooth curves and surfaces for car design [17, 18]. Bézier curves are also frequently used in computer aided graphics design [19]. The main feature of the approach is the use of control polygons where, instead of defining a curve using points on it, points near it are used. Contrary to CPI, where the curve must pass through all control points, a Bézier curve needs only be bounded by the convex hull of its control points. The theory and implementation of Bézier curve deconvolution is shown, and the concept is evaluated using simulated data, including the presence of arterial delay and dispersion and different shapes of the tissue residue function. Proof-of-concept in vivo validation is accomplished by re-evaluation of DSC-MRI data from a healthy volunteer.

Theory

Concentration and tracer kinetic theory

The assumption of a linear relationship between the change in transverse relaxation rate and concentration implies the following relationship [4, 5]:where is the MRI signal at time and is the baseline or pre-CA signal.

The AIF is the arterial concentration of contrast agent representing the intravascular tracer delivery to the local capillary network. The tissue tracer concentration is proportional to the convolution of with the tissue impulse response function, i.e., the tissue residue function R scaled by the CBF [20]:where is a constant accounting for the difference in hematocrit levels between capillaries and large vessels and the brain tissue density [21]. Estimation of the product via deconvolution provides an estimate of CBF.

MTT can be estimated as the area under the residue function [22]:and CBV is given by Eq. 4 [4]:

Bolus delay and dispersion

With a delay of between the measured AIF and the true AIF at the local tissue inlet, the measured AIF is given by , i.e., the AIF shape is retained [15]. Bolus dispersion, on the other hand, modifies the shape and amplitude of the AIF and can be described as a convolution with a vascular transport function (VTF) [10, 11]. When incorporating the effects of bolus delay and dispersion, Eq. 2 can be modified as follows:

Bézier curve deconvolution

Consider the Bernstein polynomials of degree [23]:where and . A Bézier curve of the order is a linear combination of these polynomials:where is control point . Expansion for the cubic case () yields

The control points are effectively weights that determine the influence of each basis polynomial. Equivalently, each point on the Bézier curve is a weighted average of the control points. According to Eq. 8, the Bézier curve invariably begins at and ends at . Intermediate control points need not be contained by the curve, but the curve lies wholly within the convex hull of the control points. Figure 1 provides an illustration of this concept using cubic Bézier curves with two different distributions of control points. The clear difference between the two shapes in Fig. 1a and b highlights the versatility of these curves. Figure 1c shows the third-order Bernstein polynomials that form the basis for the cubic Bézier curves. It is important to note that a Bézier curve is fully specified by its control points.

Fig. 1

Illustration of the concept of Bézier curves. a and b show cubic Bézier curves with different control point configurations and c shows the corresponding basis polynomials

In the context of deconvolution, the Bézier curve represents the residue function, and the control points constitute parameters to be determined from input data as illustrated by Fig. 2. Parameter fitting is accomplished by comparison of the modeled CTC, as depicted in Fig. 2, with the observed CTC obtained from the measured DSC-MRI signal [11]. BzD assumes a model-based approach to CBF quantification but allows a large distribution of residue function shapes within the constraints specified in the optimisation procedure.

Fig. 2

Illustration of Bézier curve deconvolution as a non-linear single-input single-output system. The residue function obtained as a cubic Bézier curve specified by the supplied control points is convolved with the measured AIF (following delay/dispersion) to give a predicted concentration time curve

Materials and methods

Implementation

Deconvolution methods were implemented in MATLAB using in-house written algorithms. For purposes of comparison, block-circulant SVD (oSVD) deconvolution was implemented according to Wu et al. [15].

BzD was implemented in a Bayesian framework adopting the maximum a posteriori (MAP) approach. The measured signal was converted to concentration using Eq. 1. Following the approach in [9], it was assumed that measurements of concentration over time followed the non-linear observation model:where is the observed concentration at time , represents the right-hand side of Eq. 2, denotes model parameters and is an observational error belonging to a zero-mean Gaussian distribution. Instead of evaluating the entire posterior distribution, optimisation was performed by seeking the mode. Particularly, the optimal parameters were determined by minimizing the negative logarithm of the factorized posterior [24].

Cubic Bézier curves were used throughout this work. Higher-order curves were associated with elevated computational expense without apparent improvement of outcome (results not shown). The residue function was forced to obey by fixing the first control point at . In addition, to ensure bounded-input, bounded-output stability, the magnitude of the last control point was fixed at zero (). The x-coordinate remained free. A penalty term was incorporated into the posterior to permit the imposing of constraints. The requirements and were enforced to obtain monotonically decreasing, bijective and positive-valued residue functions. This resulted in 5 optimisation parameters for the generation of the residue function in the form of a cubic Bézier curve as given by Eq. 8. was initially evaluated with a 10-ms sampling interval, then resampled to match the arterial input data (i.e., AIF) using cubic spline interpolation. Convolution with the measured AIF and scaling with was performed to generate a predicted CTC according to Eq. 2.

BzD with delay and dispersion modeling

Delay correction was implemented by shifting the AIF by prior to convolution with the residue function (Eq. 5). For dispersion modeling, the gamma dispersion kernel was chosen as the VTF, as suggested by Mehndiratta et al. [11]:where denotes the gamma function, determines the sharpness of the kernel and is its time-to-peak. In the limit (), the kernel approaches the Dirac delta function with no dispersion. Large values of and small values of correspond to a higher degree of dispersion [11]. The VTF was convolved with the (shifted) AIF before convolution with the residue function (Eq. 5). Arterial delay () and dispersion () were additional parameters to be determined in the optimisation.

Prior knowledge

Bayesian analysis allows the incorporation of prior knowledge into parameter estimation, which can aid in preventing physiologically implausible solutions. Gaussian distributions were chosen for each prior, with mean and standard deviation provided in Table 1. The control points were estimated on a linear scale, with a prior that permitted a wide variety of residue function shapes. A non-informative prior was chosen for CBF, which means setting a standard deviation much larger than the prior mean [11]. The prior for delay was proportional to the difference in time-to-peak between the AIF and measured CTC. Priors for the VTF were adopted from Mehndiratta et al. [11]. The algorithm was initialized at control points that favored a shape intermediate between an exponential and a boxcar residue function, and CBF was set to an arbitrary initial value of 1. Delay and dispersion parameters were initialized at their prior means.

Table 1

Prior mean and standard deviations for Bézier curve deconvolution

Parameter	Mean	Standard deviation
[P_1,x,P_1,y,P_2,x,P_2,y,P_3,x]	[8,0.5,2,0.2,15]	[8,1,4,1,100]
CBF [ml/g/s]	0.01	10^6
Delay δ [s]	-	5
VTF (s,p)	[ln2,ln2]	[2,2]

Simulations

AIF, R(t), S(t) and C(t)

A gamma-variate function was used to simulate an AIF, in accordance with previous DSC-MRI simulation studies [7-11]:where is the bolus arrival time, and . The AIF was generated over a time interval of 200 s to prevent signal truncation for the longest MTT simulated in this work (24 s). The sampling interval was .

Part of the purpose of this work was to evaluate the robustness of BzD to different underlying residue function shapes. Hence, residue functions belonging to the following family of gamma distributions were chosen [8, 25]:where with MTT = CBV/CBF and being a shape parameter. The case corresponds to an exponential residue function, = was used for the boxcar case and represented a sigmoid shape (intermediate between the exponential and the boxcar) [25]. The three different shapes are depicted in Fig. 3 for MTT = 12 s.

Fig. 3

Residue functions belonging to the gamma distribution family, for MTT = 12 s. Values of of 1, 5 and 100 represent the exponential, sigmoid and boxcar shapes, respectively

Simulations were performed for a CBV of either 4% (~ 4 ml/100 g) or 2% (~ 2 ml/100 g) representing normal gray and white matter, respectively [7, 10]. From the CTC formed by the convolution of the AIF with the residue function, signal curves were generated using Eq. 1:where was set to , and the constant was adjusted to produce a peak signal drop of at a CBF of 60 ml/100 g/min and a CBV of , resembling typical values in normal human gray matter [7, 8, 11]. The simulated AIF was also converted to signal using Eq. 13 but choosing a value of that resulted in a peak signal drop of . Zero-mean Rician noise was added to the signal curves, and all simulations were carried out using oSVD and BzD at baseline SNRs of 20 and 100 with in Eq. 2 set to 1. The oscillation index (OI) for oSVD was fixed at 0.065 for SNR = 100 and 0.035 for SNR = 20 [8]. Examples of the simulated signal data can be found in Supplementary Material 1.

No delay, no dispersion

A set of simulations was carried out without any delay or dispersion in the simulated AIF, and with the three residue function shapes shown in Fig. 3. For CBV , CBF was varied from 10 to 70 ml/100 g/min in 10 ml/100 g/min increments, and for CBV , CBF values were simulated in the range ml/100 g/min in 5 ml/100 g/min steps. As a result, the same MTT range () was obtained for both values of CBV. For each combination of , CBV, CBF and SNR, a total of 1024 CTCs were generated. The remainder of this study considered only the exponential residue function shape; this was motivated by the fact that the exponential exhibits the challenging feature of lacking an initial plateau while still being a standard shape in previous simulation studies.

Delay without dispersion

Signal curves were generated after shifting the AIF by and s relative to the simulated residue function, as in [10]. Since these shifts were not integer multiples of , the AIF was evaluated at a sampling rate of 100 ms, shifted and then resampled to . This routine was performed for CBV and CBF in the range 10–70 ml/100 g/min in 10 ml/100 g/min increments. For each combination of , CBF and SNR, 100 CTCs were formed. Deconvolution was performed using oSVD and BzD with delay correction.

Dispersion without delay

Modeling of dispersion in the absence of delay was achieved by convolving the simulated AIF with a vascular transport function. To avoid using the same VTF as the one in the deconvolution model, an exponential dispersion kernel was selected:where is a time constant for the kernel [11]. Large values of are associated with a higher degree of dispersion and in the limit , the VTF approaches the Dirac delta function with zero dispersion. Low, medium and high dispersion levels were modeled with and , respectively [10, 11]. Simulations included CBV and CBF ml/100 g/min. For each combination of dispersion level, CBF and SNR, 100 CTCs were generated, with deconvolution being performed using oSVD and BzD with dispersion correction.

Dispersion with delay

The simulated AIF was delayed then dispersed as described in the two previous sections. Low, medium and high delay and dispersion were simulated using the combination of parameters . Simulations were again performed for the same CBV, CBF range and SNR as in the previous two sections, generating 100 CTCs for each combination of these variables. Deconvolution by oSVD and BzD with correction for delay and dispersion was employed.

Determination of MTT and rMTT ratios

For each of the three AIF distortions described above (delay, dispersion and delay with dispersion), the quality of MTT estimates was evaluated using two metrics, i.e., the MTT ratio and the rMTT ratio. The MTT ratio was calculated by dividing the estimated MTT by the ground truth value. To mimic a clinical setting, where relative MTT quantification is carried out by normalizing to the MTT in a contralateral region of interest (where the AIF is assumed to be non-distorted), the rMTT ratio was introduced [11]. The rMTT ratio was, thus, defined as the estimated MTT in the presence of a distortion divided by the estimated MTT in the absence of that distortion. For example, let be the estimated MTT at a simulated AIF delay of zero and let be the MTT estimate at an AIF delay of . The rMTT ratio at the delay of is then given by .

Accuracy of delay and dispersion parameter estimation

To investigate the performance of the BzD algorithm in the estimation of the parameters describing delay and dispersion, signal curves were generated after dispersing the AIF with the same VTF as used for dispersion correction (gamma dispersion kernel) at an SNR of 100, CBV 4% and CBF 60 ml/100 g/min. Simulations were done at low, medium and high delay and dispersion, with 100 CTCs being created in each case. Different levels of dispersion were modeled by setting the gamma dispersion kernel parameters to , and for low, medium and high dispersion, respectively [11]. Analysis was carried out using BzD as described above. Relative errors in the parameters (delay), (time-to-peak of the VTF) and the sum were computed as in [11]. Note that both and were free in the fitting, but only was selected for accuracy evaluation. The results of this analysis are presented in Supplementary Materials 4 and 5.

In vivo data

A proof-of-concept evaluation of BzD was performed through a re-analysis of previously acquired DSC-MRI data from one healthy volunteer [26]. Image acquisition parameters for the DSC-MRI experiment are given in [26]. In brief, the healthy volunteer was scanned at 3 T and a global AIF was sampled from the middle cerebral artery branches in the Sylvian fissure region. A prebolus administration of a low dose of contrast agent (CA) was used to retrieve a venous output function (VOF) from the sagittal sinus prior to the main DSC-MRI experiment. After correcting the VOF area-under-curve (AUC) for the lower CA dose, a subject-specific correction factor for partial volume effects was obtained as the ratio AUC(AIF)/AUC(VOF). The same correction factor as was used by Knutsson et al. [26] was applied in the present study. The OI for oSVD was set to 0.095 [8] and in Eq. 2 was set to 0.705 [21]. Whole-brain CBF and MTT means were computed after exclusion of pixels with values above 2.5 times the mean over the entire volume as these values were assumed to originate from large vessels.

Results

Figure 4 summarizes estimated CBF against true CBF for the three different underlying residue function shapes. Each point in the plots represents the mean and standard deviation from 1024 noise realizations. BzD estimates were consistently closer to the true CBF than the corresponding oSVD estimates for the exponential residue function (). For and , BzD tended to slightly underestimate high CBF values when SNR = 20, for both CBV levels. The BzD method exhibited excellent agreement with true CBF values at SNR = 100. Underestimation of high CBF values by BzD was not evident in the boxcar () case, with even some degree of overestimation being observed at SNR = 100. While lower CBF values were generally correctly estimated by both methods, oSVD consistently failed to reproduce high CBF components (i.e., short MTT). Generally, oSVD estimates showed lower variation about their means at low SNR than the corresponding estimates with BzD. Supplementary Material 2 shows the corresponding CBF estimates, for CBV = 4% and SNR = 20, in the presence of delay and dispersion for an underlying exponential residue function.

Fig. 4

Estimated CBF plotted against true CBF for CBV = 4% and 2%, SNR = 20 and 100 and and 100, modeling the exponential, sigmoid and boxcar residue functions shapes, respectively. Results are shown for Bézier curve deconvolution and oSVD

To obtain a CBF-independent metric to conveniently compare BzD and oSVD, the CBF ratio was computed over the simulated CBF range for a given , CBV and SNR. These results are presented in Table 2, indicating that BzD outperformed oSVD given an exponential residue function and, generally, also given a sigmoid shape. For the boxcar residue function at SNR = 20, oSVD mean ratios were closer to unity than those of BzD. The performances of the two methods for boxcar shape and SNR = 100 were comparable.

Table 2

Mean CBF ratios obtained with BzD and oSVD for three different residue function shapes, SNR = 20 and 100

	SNR = 20, CBV = 4%			SNR = 100, CBV = 4%
	λ=1	λ=5	λ=100	λ=1	λ=5	λ=100
BzD	1.01±0.12	1.05±0.27	1.13±0.27	1.02±0.04	1.06±0.05	1.14±0.09
oSVD	0.69±0.16	0.89±0.20	0.99±0.24	0.80±0.12	1.02±0.11	1.13±0.09

The results show the mean and standard deviation of the CBF ratios for the ranges 10 – 70 ml/100 g/min (CBV = 4%) and 5 – 35 ml/100 g/min (CBV = 2%)

Examples of residue functions obtained with BzD and oSVD are presented in Fig. 5, for SNR = 20, MTT = 4 s (CBV = 4%, CBF = 60 ml/100 g/min) and MTT = 12 s (CBV = 4%, CBF = 20 ml/100 g/min). BzD produced smooth, monotonically decreasing, non-negative residue functions that were free of the oscillatory behavior that characterized the corresponding oSVD solutions. Moreover, BzD was in all cases in superior agreement with the true underlying residue function compared to oSVD. This latter finding is quantified more systematically in Table 3, using the root-mean-square error (RMSE) as a measure of accuracy of residue function estimation. According to Table 3, BzD provided more accurate estimates of residue function shape than oSVD in all simulated scenarios. The accuracy of residue function shape estimation with BzD declined as the underlying shape deviated from an exponential (that is, as increased); this result is coherent with the trend observed for the CBF estimation. Furthermore, for a given simulated residue function shape, the accuracy was invariably higher for higher CBV (Table 3) and for shorter MTT (Fig. 6).

Fig. 5

Typical residue functions obtained in simulation with BzD and oSVD, using an SNR of 20, and MTTs of 4 s and 12 s, for three different underlying residue function shapes

Table 3

Root-mean-square errors in residue function shape estimation, showing the mean and standard deviation over the CBF range 10–70 ml/100 g/min (CBV = 4%) and 5–35 ml/100 g/min (CBV = 2%)

	SNR = 20, CBV = 4%			SNR = 100, CBV = 4%
	λ=1	λ=5	λ=100	λ=1	λ=5	λ=100
BzD	0.03±0.02	0.04±0.05	0.06±0.06	0.02±0.01	0.02±0.02	0.04±0.04
oSVD	0.13±0.03	0.12±0.04	0.13±0.04	0.08±0.02	0.08±0.02	0.08±0.02

The results are presented for three different underlying residue function shapes, SNR = 20 and SNR = 100 with both oSVD and BzD

Fig. 6

Mean residue function fit estimated with BzD, for MTT = 4 s (CBF 60 ml/100 g/min for a CBV of 4%) and 12 s (CBF of 20 ml/100 g/min for a CBV of 4%), for SNR = 20 and 100. The dotted lines represent one standard deviation

A visual appreciation of the results in Table 3 is given in Fig. 6, which shows the mean residue function fit obtained with BzD, plus/minus one standard deviation. Results are shown for all three simulated shapes, MTT = 4 s (CBV 4%, CBF 60 ml/100 g/min) and MTT = 12 s (CBV 4%, CBF 20 ml/100 g/min), both for SNR = 20 and 100.

The accuracy of MTT estimation with BzD, BzD with delay correction and oSVD is summarized in Fig. 7. MTT and rMTT ratios were averaged over the range 3.43–24 s corresponding to CBF in the range 10–70 ml/100 g/min for a CBV of 4% and SNR = 20. The oSVD method overestimated MTT but demonstrated consistency over all simulated shifts. BzD without delay correction overestimated MTT, with its performance diminishing as the delay increased. Excellent agreement with true MTT values was achieved with BzD including delay correction. The rMTT ratios for oSVD and BzD with delay correction were close to unity across all simulated delays, consistent with the observed trend in MTT ratios. A slight underestimation was seen in the rMTT ratios for oSVD. BzD without delay correction displayed a sensitivity to delay, producing successively higher rMTT ratios as the delay increased.

Fig. 7

MTT ratio (left) and rMTT ratio (right) for oSVD, BzD and BzD with delay correction. Each bar represents the mean and standard deviation over MTT estimates for underlying CBF values in the range [10–70] ml/100 g/min, a CBV of 4% and SNR 20

The corresponding results for the case of dispersion without delay are presented in Fig. 8. MTT overestimation is apparent with oSVD, and the method exhibited a gradual decline in performance from low to high dispersion. BzD without dispersion correction displayed marked sensitivity to dispersion, showing a performance matching that of oSVD at all levels of the effect. Incorporation of a VTF into BzD drastically improved MTT estimation for low, medium and high dispersion, but led, on the other hand, to underestimation of MTT in the absence of this effect. The rMTT ratios for both oSVD and BzD with dispersion correction drifted further from unity with increasing level of dispersion. The performance of oSVD for relative MTT quantification (as captured by the rMTT ratios) was higher than that of BzD with dispersion correction, while BzD without correction for dispersion returned the worst results.

Fig. 8

MTT ratio (left) and rMTT ratio (right) for oSVD, BzD and BzD with dispersion correction. The bars show the mean and standard deviation over MTT estimates for underlying CBF values in the range [10–70] ml/100 g/min, a CBV of 4% and SNR 20

Delay and dispersion

Shown in Fig. 9 are the MTT and rMTT ratios for the case of delay coupled with dispersion. A rather similar trend to the one observed in Fig. 8 (dispersion without delay) is evident. BzD without correction for delay or dispersion exhibited lower performance when delay was combined with dispersion than with dispersion alone. When equipped with correction for both delay and dispersion, BzD appeared to outperform oSVD for absolute estimation and parallel oSVD for relative estimation. Finally, note that results similar to those shown in Figs. 7, 8, and 9 can be found in Supplementary Material 3 for an SNR of 100.

Fig. 9

MTT ratio (left) and rMTT ratio (right) for oSVD, BzD and BzD with correction for both delay and dispersion. The results show the mean and standard deviation over MTT estimates for underlying CBF values in the range [10–70] ml/100 g/min, a CBV of 4% and SNR 20

Residue functions with delay and dispersion

Evaluation of the ability to reproduce the residue function shape in the presence of delay and dispersion is summarized in Fig. 10. The plots show mean estimated residue functions for BzD (with and without correction) and oSVD, for the cases of high delay without dispersion, high dispersion without delay and high delay and dispersion. In each case, CBV was 4%, CBF 60 ml/100 g/min (typical gray matter) and SNR 20. The estimated residue function in the presence of delay and/or dispersion was characterized by a slower decay relative to the simulated shape. BzD with correction for delay and/or dispersion brought the resulting residue function into close agreement with the underlying true shape.

Fig. 10

Mean residue function fit produced with oSVD, BzD and BzD with correction for delay and/or dispersion. The plots were generated in the presence of high delay without dispersion (left), high dispersion without delay (middle) and high delay combined with high dispersion (right), at CBV = 4%, CBF = 60 ml/100 g/min and SNR = 20

Parametric maps from the analysis of clinical data with oSVD, BzD, BzD with delay correction and BzD with correction for delay and dispersion are shown in Fig. 11. BzD in any of the presented variants produced higher flow estimates than oSVD. The whole-brain CBF means for the four methods were: oSVD (42 27), BzD (58 41), BzD + delay (100 70) and BzD + delay + VTF (150 100) ml/100 g/min. Delay correction with BzD led to an increase in the calculated CBF and the incorporation of a VTF achieved the same effect.

Fig. 11

CBF, MTT and delay maps from the analysis of in vivo DSC-MRI data from a healthy volunteer using oSVD, BzD, BzD with delay correction and BzD with correction for both delay and dispersion

The MTT maps also showed variation among the four presented scenarios, with highest values being observed with oSVD (4.9 2.8) and BzD (5.3 5.8) s. MTT values for BzD methods decreased in the following order: BzD, BzD with delay correction (3.2 4.6) s, BzD with delay and dispersion correction (1.8 4.2) s. There was some discrepancy between oSVD and BzD in the delay estimation, with oSVD producing notably higher delay values than both variants of BzD shown. Interestingly, BzD with delay correction estimated higher values of delay than the same method incorporating a VTF.

Residue function shape estimation in vivo is shown in Fig. 12, for the four methods of analysis displayed in Fig. 11. The residue functions plotted are means over 4 × 4 pixels in two different ROIs. Solutions from oSVD are oscillatory, while BzD produces smooth and monotonically decaying residue functions. For both ROIs, correction for either delay or delay and dispersion results in the same effect. The residue function plots are strongly reminiscent of the results presented in Fig. 10.

Fig. 12

In vivo residue function shapes obtained with oSVD, BzD and BzD including correction for dispersion and/or delay. The plots reflect mean solutions over the 4 × 4 pixel ROIs

Discussion

Retrieval of the tissue impulse response function from the measured concentration curve and the AIF has been a persistent challenge in DSC-MRI, despite the introduction of a number of relevant mathematical methods. Recent studies include stable spline deconvolution [27], deconvolution with dispersion-compliant bases [28] assessment of the robustness of spatio-temporal deconvolution algorithms [29] as well as the introduction of machine learning approaches to DSC-MRI perfusion imaging [30]. In the present study, a non-parametric deconvolution method using Bézier curves for perfusion quantification in DSC-MRI was presented and evaluated.

In the absence of delay and dispersion, Bézier curve deconvolution showed obvious robustness to different underlying residue function shapes. The method produced accurate perfusion estimates and outperformed oSVD, except in the most extreme case (boxcar residue function at an SNR of 20). CBF estimates from BzD were accompanied by smooth, physiologically plausible tissue residue functions that obeyed the fundamental requirements of non-negative values and monotonical decrease, and they were, thus, free of the oscillations typically seen for SVD-based deconvolution methods. It is not surprising that BzD exhibited lower performance with a boxcar residue function than with the sigmoid or exponential shapes. Complete representation of such a shape (exhibiting a non-smooth behavior) using Bézier curves requires a relatively complex configuration of control points, rendering it an unlikely solution when using cubic curves. Bézier curves are inherently smooth, and thus, they favor smooth functions such as the exponential/sigmoid. It is in principle possible to bring the calculated residue function into good agreement with the boxcar shape using higher-order curves (more control points), but as highlighted earlier, this was not within the scope of the current implementation. Implementations with fourth order curves (results not shown) significantly increased computation times, and higher-order curves also had adverse effects on the accuracy of CBF estimation when combined with correction for delay and/or dispersion with an exponential underlying residue function. Incorporation of delay and dispersion correction was judged to be crucial, while accurate representation of a non-smooth residue curve had lower priority, given that non-smooth behavior is not physiologically plausible and functions resembling the exponential/sigmoid shapes are a more reasonable representation of true physiology. However, the estimates to the boxcar residue function obtained with BzD in this study appear, from visual inspection, to show closer agreement with the true shape than those obtained with the CPI method by Mehndiratta et al. [7].

BzD including a correction for delay was able to correctly reproduce MTT values in the range 3.4–24 s, for shifts of up to 6 s. Accurate estimates were also obtained by BzD with corrections for dispersion as well as for delay with dispersion, except for some degree of underestimation for zero and low levels of each effect. Zero dispersion is only achieved when the gamma dispersion kernel approaches a delta function (i.e., the sharpness parameter in the VTF must approach infinity for a time-to-peak of zero). Convergence to these values of the VTF parameters did not occur, giving rise to the observed MTT underestimation. Relative quantification results showed comparable performance between oSVD and BzD with delay and/or dispersion correction. Furthermore, oSVD overestimated MTT at all simulated shifts and dispersion levels but demonstrated a consistency that emanates from its inherent insensitivity to bolus delay. Hence, in the presence of delay and dispersion, for absolute perfusion quantification, oSVD parallels BzD when neither delay nor dispersion is corrected for. Under the same circumstances, BzD with a correction largely outperforms oSVD. For relative quantification, BzD with correction and oSVD are equally suitable, both surpassing BzD without delay or dispersion correction.

BzD without corrections provided smooth residue functions also in the presence of delay and dispersion, but both effects being signified by a spread-out solution (slower decay of the estimated residue function). Correction for delay and dispersion restored the estimated residue function to good agreement with the true shape. This result illustrates that accurate perfusion estimation with BzD, in the presence of delay and dispersion, is accompanied by physiologically plausible residue functions.

CBF maps attained from in vivo data were in general correspondence with the trends seen in the simulations, with the highest CBF values being estimated with BzD accounting for both delay and dispersion effects. Whole-brain mean CBF obtained with oSVD and BzD were in harmony with values estimated in a previous analysis of the same data [26] and with values measured with positron emission tomography [31]. However, BzD produced higher blood flow estimates with larger variability than oSVD, a result which mirrors the findings presented in Table 2. The higher spread in CBF estimates with BzD is believed to originate from low CBV regions, and we note that a more comprehensive comparison between oSVD and BzD might include a larger range of CBV values. In future work, application of both oSVD and BzD in stroke or cancer patients would also provide a more realistic comparison. BzD with delay correction, and particularly BzD with correction for both delay and dispersion, tended to overestimate CBF. This is to be expected, and it implies that low levels of delay and dispersion prevailed (cf. Figures 7, 8, 9 and Supplementary Material 5), which seems reasonable considering that the DSC-MRI data were acquired from a healthy volunteer. Very low levels of dispersion might, as observed in the simulations (cf. Figure 9), have led to an overestimation of CBF (i.e., an underestimation of MTT) when correcting for both delay and dispersion. MTT values followed the inverse of the pattern observed in CBF maps, in agreement with the central volume theorem. The result is also reasonable when viewing MTT as the area under the residue function; simulations established that both delay and dispersion correction led to a faster decay of the residue function, which is associated with a shorter MTT.

The observation that oSVD produced higher delay estimates than BzD in vivo is likely a consequence of the oscillatory nature of the oSVD residue functions coupled with the current definition of delay as the time point at which the residue function attains its maximum value. BzD with delay correction and BzD with delay and dispersion correction exhibited a notable discrepancy in delay estimation in vivo. This finding constitutes a practical example of the discussion pertaining to the results in Table S1 (Supplementary Material 4), namely that delay and dispersion are corrected in combination rather than as individual effects.

BzD also gave smooth residue functions in vivo and the effects of delay and dispersion correction reflected those observed in the simulations. An additional highlight of the in vivo results is the shape variability between the two sampled regions. It can be inferred that BzD is able to reflect in vivo residue functions with both fast and slow decay rates.

In conclusion, Bézier curve deconvolution bears substantial promise for perfusion quantification in DSC-MRI. It is a robust, non-parametric deconvolution technique that generally provides more accurate perfusion estimates than oSVD, along with realistically smooth residue functions. The power of the algorithm lies in its employment of the intrinsically smooth Bézier curves, the small number of parameters that require optimisation, and the ability to incorporate appropriate prior knowledge into the deconvolution process. We envision Bézier curve deconvolution to be of clinical value in the investigation of patients with cerebrovascular diseases, where dispersion effects are known to corrupt oSVD-based estimates of perfusion-related parameters [7, 10]. For instance, more reliable perfusion estimates obtained with the proposed approach can enable improved delineation of the salvageable ischemic penumbra in stroke patients. Furthermore, the availability of physiologically plausible residue functions opens an important avenue of research into the diagnostic potential of the residue function shape, for example, for assessment of cerebral oxygen extraction [9].

Below is the link to the electronic supplementary material.

Supplementary file1 (TIF 10654 KB) Simulated DSC-MRI data examples for MTT = 4s (top) and MTT = 24 s (bottom). The images (labeled “Image at 27 s”) show the simulated typical gray matter and white matter pixels. The data were generated at SNR = 20

Supplementary file2 (TIF 7601 KB) CBF estimates, for CBV = 4% and SNR = 20, obtained with oSVD and BzD in the presence of delay and dispersion using an exponential residue function. These results are to be compared to Fig. 4 in the main text, which shows corresponding estimates in the absence of delay or dispersion

Supplementary file3 (TIF 31707 KB) MTT and rMTT ratios from the study of delay and dispersion at SNR = 100. Similar results are shown for an SNR of 20 in Figs 7–9 in the main text

Supplementary file4 (PDF) Accuracy of delay and dispersion parameter estimation

Supplementary file5 (TIF 23053 KB) MTT estimates from BzD and oSVD obtained from the study of delay correction in the presence of dispersion without delay. Notice how BzD with delay correction improves MTT estimates even though delay is absent