Loading...

Summer Research Fellowship Programme of India's Science Academies

Joint estimation of Arterial Input Function and tracer kinetic parameters from undersampled DCE-MRI

Aparna Karnik

3rd year Bachelor of Engineering, Electronics and Communication department, Ramaiah Institue of Technology, Bangalore 560054

Dr. Phaneendra Kumar Yalavarthy

Associate Professor, Department of Computational and Data Sciences, IISc, Bangalore 560012

Abstract

Dynamically Contrast Enhanced- Magnetic Resonance Imaging (DCE-MRI) is a relatively new modality that helps in both pre-clinical tumor modelling as well as clinical oncology. It requires acquisition of MR images before and during the administration of contrast agent (CA) or tracer. Further quantitative analysis of the acquired data provides information for diagnosis and monitoring response to medication. In DCE-MR images we need high spatial and temporal resolution. However, a tradeoff exists between the two. In order to obtain a large spatial coverage at high temporal resolution, undersampling technique is a feasible option. This brings out the need for a good reconstruction technique to make up the loss of data in the k-space. AIF and TK parameters can be jointly estimated from the undersampled k-space data directly without reconstructing the image of the tissue. This work aims at the comparison of Lp-norm regularisation schemes (p≤1) and L2-norm regularisation scheme to enforce model consistency and compare the accuracy and the reconstruction time of the two methods with 20x, 40x, 60x, 100x and 120x undersampled data.

Keywords: T1 weighted image, k-space, Arterial Input Function, tracer kinetic parameter

Abbreviations

Abbreviations
DCE-MRIDynamic Contrast Enhanced Magnetic Resonance Imaging
 CAContrast Agent
AIF  Arterial Input Function
 TK Tracer Kinetic
 ROI Region Of Interest
pat-AIF patient specific Arterial Input Function 
 pop-AIFpopulation average Arterial Input Function 

INTRODUCTION

Background

A tumor, however benign it might be, cannot depend on diffusion metabolics for its proliferation. Thus, it starts developing its own vasculature which is known as angiogenesis. These new vessels, unlike the vasculature of normal healthy tissues, are leaky, missing epithelial cells and have blind ends. DCE-MRI exploits these differences to characterise the neo-vasculature, indirectly the tumor, as they appear as a bright spot due to leakage of the CA.

angiogenesis_1.jpg
    Angiogenesis . Source: http://sites.google.com/site/bmed360

    The conventional Nyquist sampling is unable to simultaneously provide the adequate temporal resolution and spatial coverage. DCE-MR images acquired at Nyquist rate can provide around 5-10 slices in the time span of flow of CA. To achieve the required temporal resolution at Nyquist rate, spatial resolution can be coarsened, but this results in a compromise with the quality of the image and also hinders the clinical evaluation of the image. Thus techniques involving undersampling and reconstruction have been developed to simultaneously obtain the required spatial coverage and temporal resolution. A good reconstruction model uses the prior information provided by standard modules on the tracer kinetics in capillaries while complying to the model consistency constraints and data fidelity. Early work on reconstruction was based on reconstructing the MR image from the under-sampled k,t-space and then applying TK modelling software to obtain the high resolution TK maps based on the reconstructed images. However, newer studies have been focusing on directly using the TK models to reconstruct the undersampled k,t-space data. These studies have shown to allow higher undersampling rates.

    Statement of the Problems

    Studies have been using L2-norm regularisation to enforce model consistency. Recent studies have also suggested the possible use of Lp-norm (p≤1) regularisation methods. This review draws an efficient comparison between the two methods of regularisation and also tries to compare the existing backward model methods i.e. Patlak model and eTofts model.

    LITERATURE REVIEW

    Dynamically Contrast Enhanced (DCE) - MRI

    DCE-MRI is the serial acquisition of T1 weighted images, before and after the CA is intravenously administered, in the ROI. CAs are usually paramagnetic material which decrease T1 time. The addition of CA makes the tumors in the ROI to appear brighter than the surrounding healthy tissue. MR signal intensity of any tissue is directly dependant on the local concentration of CA. By this relation, due to angiogenesis, the concentration of CA is much higher in a malignant tumor than in a healthy tissue. Upon quantitatively assessing the signal intensity in the time span from the administration of CA to its ejection, an insight into parameters related to blood flow, vessel permeability, tissue volume fractions etc., can be obtained.

    To support the DCE-MRI data quantitative analysis, compartment models are used. Human body can be represented as one or more 'compartments' in which the CA dynamically flows. A compartment is a bound space that a CA can occupy and whose volume remains constant on the time scale of the experiment. In each compartment, it is assumed that, the CA is uniformly distributed throughout, without any time delay.

    The simplest as seen in fig 1 from ​Thomas E. Yankeelov, 2007 is the single compartment method wherein it is assumed that the entire human body is one single compartment and the drug under consideration reaches an equilibrium throughout the body. This assumption is seldom held good as CAs are typically distributed with different rates in different tissues.

    single_compartment.jpg
      The simplest pharmacokinetic model is the single compartment model depicted above. A drug enters the compartment, V, with a specific rate constant ka, and is eliminated with a rate constant ke. Source:  Thomas E. Yankeelov, 2007

      A more realistic model is the two compartment model which considers the extracellular intravascular space i.e. the blood plasma to be the central compartment and the extracellular extravascular space as the peripheral compartment.

      two_compartment.jpg
        A more realistic approach is described by the two-compartment model which considers the extracellular intravascular space (blood plasma) to be the central compartment (V1), and the extravascular-extracellular space (EES) to be the peripheral compartment (V2). In this model, the CA is introduced into the vasculature and transported into the EES in a reversible process characterized by a distribution rate constant (k12) and a redistribution rate constant (k21). Source:  Thomas E. Yankeelov, 2007

        Arterial Input Function

        Obtaining the Arterial Input Function (AIF) is by far the most technically demanding process in the acquisition of DCE-MRI. The AIF is the estimate of the concentration of CA in the blood plasma as a function of time, Cp(t). There are three main types of obtaining the AIF for data analysis.

        The first method involves introducing an arterial catheter into the subject and measuring the concentration of CA during the imaging process. This method is highly accurate but at the same time it is invasive and has very poor temporal resolution.

        Patient specific AIF can be obtained from the DCE-MRI data set itself. The signal intensity changes, due to the passing CA, in blood plasma and the surrounding tissues is simultaneously measured and later converted to intravascular concentration levels of CA using methods discussed later. The advantage of this method is that AIF is obtained at an individual basis while being completely non-invasive.

        The third method, population average AIF, assumes that the AIF is similar for similar subjects. The AIF is measured for a small cohort of subjects and the same is assumed to be the valid AIF for the later subjects also. The merit of this method is its simplicity in both acquisition and analysis. The major drawback is that the changes introduced by the pathology under study are ignored.

        Backward model

        The concentration change at every pixel can be used to measure the flow of CA at that pixel. There are various models that can be used to measure the flow of CA at a pixel. Some of the most common ones are the Patlak model and the eTofts model.

        Patlak model

        The Patlak model is a linear model which follows the two compartment theory.

        Ct(t)=P(θ)=Vp.Cp(t)+Ktrans0tCp(τ)dτC_t(t)=P(\theta)=V_p.C_p(t)+K^{trans}\int_0^tC_p(\tau)\operatorname d\tau

        where Ct is the total CA concentration, Vp is the fractional plasma volume per tissue volume, Cp is the AIF and Ktrans is the forward transfer model. Since this is a linear model, concentration can be retrieved by taking the pseudo inverse.

        θ=P1(C)\theta=P^{-1}\left(C\right)

        eTofts model

        The extended Tofts model assumes a term Kep which accounts for the backward flow of fluids into capillaries.

        Ct(t)=P(θ)=Vp.Cp(t)+Ktrans0tCp(τ)e-Kep(t-τ)dτ

        Since, eTofts model is non-linear in nature, a simple pseudo-inverse technique is not sufficient to obtain concentration levels back. An iterative algorithm can be used to solve this model fitting:

        θ=argminθP(θ)C22\theta=\underset\theta{argmin}\left\|P(\theta)-C\right\|_2^2

        As the name suggests, Forward model is essentially the reverse operation of Backward model, i.e., the concentration maps are derived from the TK parametric maps.

        Ct=P(θ,Cp)C_t=P(\theta,C_p)

        METHODOLOGY

        Concepts

        Model consistency constraint

        Guo et. al, 2016 and Megha propose jointly estimating contrast concentration versus time images (C) and TK parameter maps (θ) from undersampled data (y) by solving the least square problem:

        (C,θ)=argminC,θUFE(ψC+S0)y22+λ.P(θk)CPP(C,\theta)=\underset{C,\theta}{argmin}\left\|UFE(\psi C+S_0)-y\right\|_2^2+\lambda.\left\|P(\theta^k)-C\right\|_P^P

        The first term, L2 norm, represents data consistency, i.e., C should be consistent with the measured data y by signal equation (ψ), undersampling mask (U), fourier transform (F) and sensitivity encoding (E). S0 is the first time frame image which is fully sampled. S0 is usually obtained before administration of the CA. The second term, LP norm (​Pant et. al., 2014​) represents the model consistency constraint, where, C must be consistent with the modelling TK parameter maps i.e. Patlak or eTofts in this scenario.

        This review is on joint estimation of AIF and TK parameters. AIF can be incorporated in equation (5),

        (C,θ,AIF)=argminC,θ,AIFUFE(ψC+S0)y22+λP(θ,AIF)Cpp(C,\theta,AIF)=\underset{C,\theta,AIF}{argmin}\left\|UFE(\psi C+S_0)-y\right\|_2^2+\lambda\left\|P(\theta,AIF)-C\right\|_p^p

        we solve each variable alternatively as shown in (8), (nth iteration):

        θn+1,AIFn=P1(Cn+1)\theta^{n+1},AIF^n=P^{-1}(C^{n+1})

        equation (8) is backward TK modelling, which also estimates patient specific AIF. For the estimation of pat-AIF, the arterial ROI must be identified.

        Methods

        The aim of this project was to do a comparative analysis of the methods proposed by Guo et. al., and Megha Goel, i.e., L2 norm regularisation and Lp norm regularisation (p≤1) in model consistency constraints using both Patlak and eTofts model.

        Data-set

        The data used in this project is the In-vivo dataset, which has been made available by Guo et. al, 2016. This consists of nine fully sampled brain tumor DCE- MR images of a male patient of age 65 years diagnosed with glioblastoma. The parameters of this data set are as follows: field of view = 22 x 22x 4.2 cm3, spatial resolution = 0.9 x 1.3 x 7.0 mm3, temporal resolution = 5s, 50 time frames, eight receiver coils, flip angle = 15°, echo time = 1.3ms, repetition time = 6 ms, DEPOSIT1 has been performed before DCE-MRI, with flip angle of 2°, 5°, and 10° to estimate precontrast T1 and M0 maps. The contrast agent that has been used is gadobenate dimeglumine.

        The data was retrospectively undersampled using randomised Golden angle Ratio at R= 20x, 40x, 80x, 100x and 120x. Gaussian noise was added to the image. For patient specific AIF, the ROI was detected and AIF was calculated by applying the forward model to the fully sampled image S0. A modified version of Conjugate gradient method proposed by ​Pant et. al., 2014​​ was used to to optimize the function.

        RMSE

        Root Mean Square Error is used as the figure of merit in this review to compare the different methods and models used. The Ktrans maps obtained by reconstruction of the undersampled data is quantitatively compared with the Ktrans maps obtained from the fully sampled S0 pixel by pixel, in the ROI.

        Normalised RMSE is calculated as

        nRMSEKTrans=1Kfs,90trans.i=0nxm(KfullysampledTransiKundersampledTransi)2n    mnRMSE_{K^{Trans}}=\frac1{K_{f_s,90}^{trans}}.\sqrt{\frac{{\displaystyle\sum_{i=0}^{nxm}}\left(K_{fully-sampled}^{Trans_i}-K_{under-sampled}^{Trans_i}\right)^2}{n\;*\;m}}

        RMSE was also calculated for the AIF that was estimated in the models that used pat-AIF. The AIF that was estimated was quantitatively compared with the AIF obtained from the fully sampled S0, for every time frame.

        nRMSEAIF=1AIFfs,90i=050(AIFifully  sampledAIFiundersampled)250nRMSE_{AIF}=\frac1{AIF_{f_s,90}}\sqrt{\frac{\sum_{i=0}^{50}\left(AIF_i^{fully\;sampled}-AIF_i^{under-sampled}\right)^2}{50}}

        RESULTS AND DISCUSSION

        tumors_2.jpg
          Reconstruction of the Ktrans maps for R=20x, 40x, 60x, 100x and 120x for L1, L2 and L0.7 norm regularisation. The reconstructed maps are in the tumor ROI 

          Fig 4 shows the Ktrans maps of the ROI for L1 and L0.7 norms using Patlak model and L2 using eTofts.

          rms_values_3.jpg
            Graph depicting nRMSE values of L2 norm regularisation for Patlak model of pat-AIF and pop-AIF and eTofts  of pop-AIF, L1 norm and L0.7 norm regularisation of Patlak Pat-AIF.

            Fig 5 and table 1 show the RMSE values for the undersampling values 20x, 40x, 60x, 100x and 120x which were performed for Patlak and eTofts model using L2-norm, L1-norm and L0.7-norm regularisation on model consistecy constraint. As seen on the graph, L0.7- norm regularisation and L1-norm pat-specific have performed fairly well even at 120x undersampling when compared to the other methods.

            nRMSE values for R=20x, 40x, 60x 100x and 120x, for  L2 norm regularisation for Patlak and eTofts models using pat-AIF and pop-AIF, and Patlak model using pat-AIF with L1 and L0.7 norm regularisation.
             Model/Undersampling rate20  40 60 100 120
             L2 Patlak pat-specific0.1677 0.2041 0.2419 0.2962 0.2755
             L2 Patlak pop-avg0.362 0.3502 0.3895 0.4105 0.4114
             L1 Patlak pat-specific0.1037  0.1349 0.1553 0.2011 0.2548
             L0.7 Patlak patient specific0.1023  0.13560.1601  0.2102 0.2789
             L2 eTofts Pop-avg 0.2786 0.28810.2442  0.2482 0.2568

            eTofts has maintained a steady average RMSE of 0.26 throughout all the undersampling rates. Whereas, L2-norm regularisation of Patlak model using pop-AIF has performed poorly having an average of 0.38 RMSE worsening with progressing undersampling rates. L1 and L0.7 norm regularisations have performed exceptionally well with negligible error rates.

            AIF_1.jpg
              graphs of AIF calculated from the fully sampled MR image and the estimated AIF obtained from the data for L0.7, L1 and L2 norm regularisations for the 50 slices of images. All three graphs are for undersampling rate R=60x

              Fig 6 depicts the graph of AIF values for the 50 slices of images in the data for L0.7, L1 and L2 regularisations at R=60x undersampling. As the value of p in Lp-norm reduces, there is a significant rise in accuracy in the estimation of the AIF peak. Conversely, as p value increases, the average estimation has improved of the AIF after the 15th slice. Hence, it can be said that a trade-off exists between the two in terms of estimation of AIF.

              rms_values_AIF_1.jpg
                nRMSE values of AIF estimated for the pat-AIF using L1, L2 and L0.7 norm regularisation for undersampling values R=20x, 40x, 60x, 100x and 120x

                Fig7 shows the RMSE values of AIF estimation for different undersampling rates. As seen, L1 norm is the least erroneous of the three. Although, for this range of undersampling rates, all three methods have minuscule error rates.

                As suggested by ​Khalifa et al., 2014​, we can observe the simplicity of Patlak model as it omits reverse vascular transfer coefficient, Kep. This highly simplifies the model and makes it linear thus facilitating the easy extraction of concentration values from TK parameters. At the same time eTofts model, even though slightly more accurate than Patlak in some cases, fails to perform adequately due to its acceptance of pop-AIF and less flexibility in comparison with Patlak model. For all the above reasons, the Patlak model has evolved to be more popular in the quantitative analysis of DCE-MRI than eTofts model. Studies also claim that the eTofts method sometimes produces drastic errors as opposed to Patlak model.

                CONCLUSION

                We explored the implementation of the proposed methods of joint estimation of AIF and TK-parameters, weighing them based on their RMS errors. While L2-norm regularisation of Patlak model failed to perform on heavily undersampled data, the same on eTofts model seemed to show stable performance but with low accuracy. As we moved to Lp-norm (p≤1) , the accuracy visibly improved in TK parameter estimation, but at the cost of increased error in AIF estimation.

                Methods which strike a balance between the L2 and Lp norms incorporating the strengths of each and also reaching the optimum solution quicker are worth exploring. The scarcity of DCE-MRI datasets that provide complete AIF information is a major setback in this regard.

                Bibliography

                • Fox, Leslie.An introduction to numerical linear algebra. No. 04; QA251, F6.. Oxford: Clarendon Press, 1964.
                • Joint Estimation of AIF and TK-parameters from Under-sampled DCE-MRI k-t space Data Using Sparse Recovery Methods, MTech Thesis, Megha Goyal
                • Heggie, John CP. "Magnetic resonance imaging: Principles, methods and techniques by perry sprawls."Australasian Physical & Engineering Science in Medicine24.1 (2001): 57-57.

                References

                • Thomas Yankeelov, John Gore, 2007, Dynamic Contrast Enhanced Magnetic Resonance Imaging in Oncology:Theory, Data Acquisition,Analysis, and Examples, Current Medical Imaging Reviews, vol. 3, no. 2, pp. 91-107

                • Yi Guo, Sajan Goud Lingala, Yinghua Zhu, R. Marc Lebel, Krishna S. Nayak, 2016, Direct estimation of tracer-kinetic parameter maps from highly undersampled brain dynamic contrast enhanced MRI, Magnetic Resonance in Medicine, vol. 78, no. 4, pp. 1566-1578

                • Jeevan K. Pant, Wu-Sheng Lu, Andreas Antoniou, 2014, New Improved Algorithms for Compressive Sensing Based on $\ell_{p}$ Norm, IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 61, no. 3, pp. 198-202

                • Fahmi Khalifa, Ahmed Soliman, Ayman El-Baz, Mohamed Abou El-Ghar, Tarek El-Diasty, Georgy Gimel'farb, Rosemary Ouseph, Amy C. Dwyer, 2014, Models and methods for analyzing DCE-MRI: A review, Medical Physics, vol. 41, no. 12, pp. 124301

                Source

                • Fig 1: http://sites.google.com/site/bmed360
                • Fig 2: Yankeelov, Thomas E., and John C. Gore. "Dynamic contrast enhanced magnetic resonance imaging in oncology: theory, data acquisition, analysis, and examples." Current medical imaging reviews 3.2 (2007): 91-107.
                • Fig 3: Yankeelov, Thomas E., and John C. Gore. "Dynamic contrast enhanced magnetic resonance imaging in oncology: theory, data acquisition, analysis, and examples." Current medical imaging reviews 3.2 (2007): 91-107.
                More
                Written, reviewed, revised, proofed and published with