To propose and assess a new method for model-free quantification in perfusion MRI.

Estimation of haemodynamic properties using bolus-tracking perfusion MRI is typically accomplished by deconvolution of the tissue concentration curve ($$$C_t$$$) with an arterial input function (AIF) in order to yield the tissue impulse response function (TRF), which is the residue function scaled with tissue perfusion , i.e., $$$TRF(t)=f_t\cdot R(t)$$$ . The fact that deconvolution of noisy signals is an ill-posed and ill-conditioned problem has led to the proposal of many different deconvolution methods in perfusion MRI. The most commonly used is the singular value decomposition (SVD) based method¹, often accompanied by fix or adaptive truncation. To produce more physiologically plausible residue functions without spurious oscillations, methods based on, e.g., Tikhonov regularization², stochastic processes³, and control point interpolation⁴ (CPI), the latter two implemented in Bayesian frameworks, have also been proposed.

In this work, we introduce a deconvolution method based on Bézier curves⁵, which are inherently continuous and smooth in general. The concept of Bézier curves was developed by de Casteljau (at Citroën), and independently by Bézier (at Rénault) during the 1960s, to produce smooth curves and surfaces for car design. Bézier curves are also common in computer graphics, for example, in vector graphics software. The proposed deconvolution method, here referred to as ‘Bézier deconvolution’ (BzD), resembles the CPI method⁴, although BzD employs control polygons rather than pure control points. This means that the curve does not pass through all control points, although the curve is always contained within the convex hull of the control points. Due to their smooth characteristics and the possibility to control their shape, we believe that Bézier curves could allow for physiologically realistic residue function shapes with rather few parameters. In this work, we demonstrate initial experiences from the application of BzD to dynamic susceptibility contrast MRI (DSC-MRI) in vivo data.

A Bézier curve of the n^th order is given by

$$\mathbf{B}_n(\tau)=\sum_{i=0}^n\left(\begin{array}{c}n\\ i\end{array}\right)\left(1-\tau\right)^{n-1}\tau^i\mathbf{P}_i$$

where $$$\tau \in [0,1]$$$ and $$$\mathbf{P}_i$$$ denotes the coordinates of control point $$$i$$$. Cubic (B₃) and quartic (B₄) Bézier curves (see Fig. 1) can thus be expressed as follows:

$$\mathbf{B}_3(\tau)=(1-\tau)^3\mathbf{P}_0+3(1-\tau)^2\tau\mathbf{P}_1+3(1-\tau)\tau^2\mathbf{P}_2+\tau^3\mathbf{P}_3$$

$$\mathbf{B}_4(\tau)=(1-\tau)^4\mathbf{P}_0+4(1-\tau)^3\tau\mathbf{P}_1+6(1-\tau)^2\tau^2 \mathbf{P}_2+4(1-\tau)\tau^3\mathbf{P}_3+\tau^4\mathbf{P}_4$$

$$$\mathbf{P}_0=\left\{0,1\right\}$$$ and $$$P_{n,y}=0$$$ were fixed to ensure $$$R(0)=1$$$ as well as bounded-input, bounded-output stability (see Fig. 1c). To obtain a non-negative monotonically decreasing one-to-one solution, the conditions $$$0\leq P_{i,y}\leq 1$$$ and $$$0\leq P_{i,x}\leq P_{n,x}$$$ were also required. This resulted in 5 and 7 parameters for the cubic and quartic Bézier curves, respectively. The solution was resampled at the time points of the measured data before convolution with the AIF. The algorithm was implemented in a Bayesian framework⁶, using a non-informative prior for the CBF.

For proof of concept, the BzD method was applied to DSC-MRI data from a healthy volunteer, acquired with the following protocol: Single-shot gradient-echo EPI, matrix=128×128, voxel size=1.72×1.72×5 mm³, 20 slices, FA=60º and acquisition time=1:30 min. The project was approved by the local ethics committee, and written informed consent was obtained. A global AIF ($$$C_{AIF}$$$) was selected automatically, and the TRF was estimated based on the following model: $$$C_t(t)=\kappa f_t\left[C_{AIF}(t+\delta)\otimes R(t)\right]$$$, where $$$\kappa$$$ is a hematocrit and tissue density scaling factor (set to 0.701) and $$$\delta$$$ is a parameter taking into account the delay between tissue and arterial input curves. Mean transit time (MTT) was estimated as the time-integral of the residue function.

Figure 2 displays parametric maps from the in vivo analysis, obtained using BzD and oSVD. The CBF values were higher with BzD, compared to oSVD. MTT and delay maps were substantially different between BzD and oSVD, and oSVD produced generally longer MTT and delay estimates. Cubic and quartic BzD produced very similar results overall. Figure 3 shows density plots of all residue functions from the slice in Figure 2. BzD produced more realistic residue function shapes than oSVD.

We implemented a new deconvolution method based on Bézier curves for quantitative perfusion MRI. The BzD method produced higher perfusion values than the commonly used oSVD method (which is prone to non-physiological residue function shapes and underestimation due to truncation of singular values). The MTT and delay maps were also markedly different between BzD and oSVD.

A systematic comparison between BzD and the CPI method⁴ would be of considerable interest, to assess whether control polygons show advantages over control points. Furthermore, realistic model-free residue function shapes can also yield capillary transit time distributions, which can be used for model-free estimation of oxygen extraction capacity and metabolic rate of oxygen capacity⁷ (by application of the capillary transit time heterogeneity model⁸).