Perfusion, which represents the capillary blood flow, may be measured using dynamic susceptibility contrast MR Imaging (DSC-MRI)¹. The DSC-MR Image analysis is typically performed using time-series of T2*-weighted or T2-weighted single-shot echo-planar images (EPI). This base image series may contain several hundreds of images acquired in quick succession to provide a dynamic display (one scan every 1 to 2 seconds over 1 to 2 minutes). This tissue concentration of contrast agent obtained using this dynamic series should be deconvolved with concentration in the supplying blood vessels to obtain the response of an ideal input¹. This input is known to be residue function (RF) from which the quantitative maps, such as cerebral blood volume (CBV), cerebral blood flow (CBF), mean transit time (MTT), and time until the peak of residue function (T_max) maps are computed¹. In this work, the QR based decomposition of convolution matrix (obtained from concentration in the supplying blood vessels) combined with traditional Tikhonov regularization is proposed to accurately estimate the deconvolved response (RF).

In DSC-MRI, a mathematical description in the matrix form that is necessary to describe the complex relationship between signal intensity and contrast agent concentration by taking into account the physical contrast mechanisms, has the form $$C_t = \mathbf{C_a}R + e,$$ where $$$C_t$$$ is the vector samples tissue concentration $$$C(t_1, t_2, ..,t_n)^T$$$ and $$$R$$$ (known as the tissue residue function from which the quantitative perfusion maps are obtained) is given by $$$(R(0), …R(t_n))^T$$$, where $$$n$$$ is the number of time points and $$$T$$$ denotes the transpose operation with $$$e$$$ representing the noise vector. The entries of $$$\mathbf{C_a}$$$ (dimension: $$$n\times n$$$) are in the circulant form formed by the arterial input function. The matrix $$$\mathbf{C_a}$$$ is nearly singular (highly ill-conditioned) in nature. The e represents the additive noise in the MR signal. To obtain R one needs to minimize the following objective function¹: $$\Omega = || \mathbf{C_a}R - C_t||^2 + \lambda^2 ||R||^2.$$ Traditionally this is minimized either using Oscillatory-limited SVD (oSVD)² or Frequency-Domain deconvolution (FDD)³.

The proposed method performs QR decomposition of $$$\mathbf{C_a}$$$ by Lanczos bidiagonalization as given in Ref. [4]. With these bidiagonalizations the problem of deconvolution decomposes into the iterative Least Squares QR. It turns $$\mathbf{C_a}R - C_t = U_{k+1}(B_kr_k - β_0e_1); R = \mathbf{V_k}r_k.$$ Here $$$\mathbf{B}$$$ represents the lower bidiagonal matrix (dimension $$$k+1\times k$$$), with $$$(\alpha_1,…,\alpha_k)$$$ on the main diagonal and $$$(\beta_1,…, \beta_k)$$$ in the lower sub-diagonal), $$$\beta_0$$$ is the $$$L_2$$$-norm of the $$$C_t$$$, $$$\mathbf{U}$$$ and $$$\mathbf{V}$$$ represent the left and right orthogonal Lanczos matrices, respectively. The unit vector of dimension $$$k\times 1$$$ is represented by $$$e_k$$$ (=1 at the k^th row and 0 elsewhere). The dimensions of $$$\mathbf{U_k}$$$ and $$$\mathbf{V_k}$$$ are $$$(n\times k)$$$ and $$$(n\times k)$$$, with $$$k$$$ representing the number of iterations the bidiagonalization is performed. Equivalently, it converts the minimization function into $$\Omega = || \mathbf{B_k}r_k - \beta_0e_1||^2 + \lambda^2 ||r_k||^2.$$ This results in update equation $$r_k = (\mathbf{B_k}^T\mathbf{B_k} + \lambda^2\mathbf{I})^{-1} β_0\mathbf{B_k}^Te_1); R = \mathbf{V_k}r_k.$$ Note that $$$k \ll n$$$ leading to dimensionality reduction of the original problem (i.e. lesser number of linear system of equations to be solved)⁴.

The $$$k$$$ is chosen based on the oscillation index of the deconvolved signal similar to oSVD approach², giving it the name Oscillatory-limited QR method (oQR).

The open-source digital phantom data was used for evaluating the performance of discussed deconvolution methods, including the proposed one. The digital phantom dataset contained 16 slices, with slice-1 containing the AIF and VOF⁵. Three types of residue functions, exponential, linear, and box-shaped, were simulated spanning slices 2-6, 7-11, and 12-16 respectively. The forms of these residue functions are given in the Appendix-E1 of Ref. [5]. The SNR of the data is maintained at 40. The computed root mean square error for perfusion maps and expected perfusion using the discussed methods are given in the Table-1. Clearly, in the MTT for exponential and linear type, the proposed method (oQR) yields better results. For exponential type residue function, the oQR results are superior compared to other methods. For Box type residue function, the FDD provides best performance with oQR faring better than oSVD. The same comparison using the stroke patient with diagnosis of right Middle Cerebral Artery (MCA) stroke DSC-MRI (1.5 Tesla) data is given in the Fig.1. Even here, the stroke region, core along with penumbra, indicated with red arrow yields superior results with the proposed method (oQR). In summary, a superior and novel deconvolution approach based on QR decomposition in Tikhonov regularization framework has been developed for post-processing of DSC-MRI data.