Direct Estimation of Pharmacokinetic Parameters from DCEMRI using Deep CNN with Forward Physical Model Loss
Abstract
Dynamic contrastenhanced (DCE) MRI is an evolving imaging technique that provides a quantitative measure of pharmacokinetic (PK) parameters in body tissues, in which series of weighted images are collected following the administration of a paramagnetic contrast agent. Unfortunately, in many applications, conventional clinical DCEMRI suffers from low spatiotemporal resolution and insufficient volume coverage. In this paper, we propose a novel deep learning based approach to directly estimate the PK parameters from undersampled DCEMRI data. Specifically, we design a custom loss function where we incorporate a forward physical model that relates the PK parameters to corrupted imagetime series obtained due to subsampling in kspace. This allows the network to directly exploit the knowledge of true contrast agent kinetics in the training phase, and hence provide more accurate restoration of PK parameters. Experiments on clinical brain DCE datasets demonstrate the efficacy of our approach in terms of fidelity of PK parameter reconstruction and significantly faster parameter inference compared to a modelbased iterative reconstruction method.
1 Introduction
Dynamic contrastenhanced (DCE) MRI involves the administration of a shortening Gadoliniumbased contrast agent (CA), followed by the acquisition of successive weighted images as the contrast bolus enters and subsequently leaves the organ [9]. In DCEMRI, changes in CA concentration are derived from changes in signal intensity over time, then regressed to estimate pharmacokinetic (PK) parameters related to vascular permeability and tissue perfusion [6]. Since perfusion and permeability are typically affected in the presence of vascular and cellular irregularities, DCE imaging has been considered as a promising tool for clinical diagnostics of brain tumours, multiple sclerosis lesions, and neurological disorders where disruption of bloodbrain barrier (BBB) occurs. [7, 4].
Despite its effectiveness in quantitative assessment of microvascular properties, conventional DCEMRI is challenged by suboptimal image acquisition that severely restricts the spatiotemporal resolution and volume coverage [2, 3]. The shortest possible scanning time often leads to limited spatial resolution hampering detection of small image features and accurate tumor boundaries. Low temporal resolution hinders accurate fitting of PK parameters. Furthermore, volume coverage is usually inadequate to cover the known pathology, for instance in the case multiple metastatic lesions [3]. Facing such severe constraints, DCE imaging can significantly benefit from undersampled acquisitions.
So far, existing works in [6, 2, 11] have proposed compressed sensing and parallel imaging based reconstruction schemes to accelerate DCEMRI acquisitions, mainly targeting to achieve better spatial resolution and volume coverage while retaining the same temporal resolution. These methods are referred to as indirect methods [3] because they are based on the reconstruction of dynamic DCE image series first, followed by a separate step for fitting the PK parameters on a voxelbyvoxel level using a tracer kinetic model [9]. More recently, a modelbased direct reconstruction model [3] has been proposed to directly estimate PK parameters from undersampled (k,t) space data. The direct reconstruction method generally poses the estimation of PK maps as an error minimization problem. This approach has been shown to produce superior PK parameter maps and allows for higher acceleration compared to indirect methods. However, the main drawback of this method is that parameter reconstruction of an entire volume requires considerably high computation time.
Motivated by the recent advances of deep learning in medical imaging, in this paper, we present a novel deep learning based approach to directly estimate PK parameters from undersampled DCEMRI data. First, our proposed network takes the corrupted imagetime series as input and residual parameter maps, which represent deviations from a kinetic model fitting on fullysampled imagetime series, as output, and aims at learning a nonlinear mapping between them. Our motivation for learning the residual PK maps is based on the observation that residual maps are more sparse and topologically less complex compared to target parameter maps. Second, we propose the forward physical model loss, a custom loss function in which we exploit the physical relation between true contrast agent kinetics and measured timeresolved DCE signals when training our network. Third, we validate our method experimentally on human in vivo brain DCEMRI dataset. We demonstrate the superior performance of our method in terms of parameter reconstruction accuracy and significantly faster estimation of parameters during testing, taking approximately 1.5 seconds on an entire 3D test volume. To the best of our knowledge, we present the first work leveraging the machine learning algorithms – specifically deep learning – to directly estimate PK parameters from undersampled DCEMRI timeseries.
2 Methods
We treat the parameter inference from undersampled data in DCE imaging as a mapping problem between the corrupted intensitytime series and residual parameter maps where the underlying mapping is learned using deep convolutional neural networks (CNNs). We provide a summary of general tracer kinetic models applied in DCEMRI in Sec. 2.1, formulate the forward physical model relating the PK parameters to undersampled data in Sec. 2.2, finally describe our proposed deep learning methodology for PK parameter inference in Sec. 2.3.
2.1 Tracer Kinetic Modeling in DCEMRI
Tracer kinetic modeling aims at providing a link between the tissue signal enhancement and the physiological or socalled pharmacokinetic parameters, including the fractional plasma volume (), the fractional interstitial volume (), and the volume transfer rate () at which contrast agent (CA) is delivered to the extravascular extracellular space (EES).
One of the wellestablished tracer kinetic models is known as Patlak model [8]. This model describes a highly perfused two compartment tissue, ignoring backflux from the EES into the blood plasma compartment. The CA concentration in the tissues is determined by,
(1) 
where represent image domain spatial coordinates, is the CA concentration over time, and denotes the arterial input function (AIF) which is usually measured from voxels in a feeding artery.
In this work, we specifically employ the Patlak model for tracer pharmacokinetic modeling and estimation of ground truth tissue parameters. This model is a perfect match for our DCE dataset because it is often applied when the temporal resolution is too low to measure the cerebral blood flow, and it has been commonly used to measure the BBB leakage with DCEMRI in acute brain stroke and dementia [4, 9]. An attractive feature of Patlak model is that the model equation in (1) can be linearized and fitted using linear least squares which has a closedform solution, hence parameter estimation is fast [9].
2.2 Forward Physical Model: From PK Parameters to Undersampled Data
Figure 1 depicts the conventional and forward model approaches relating the PK parameter estimation to undersampled or fullysampled kspace data, and vice versa. For direct estimation of PK parameters from the measured kspace data, as proposed in [1, 3], a forward model can be formulated by inverting the steps in the conventional model as follows:

Given the sets of PK parameter pairs () and arterial input function , CA concentration curves over time are estimated using the Patlak model equation in (1).

Dynamic DCE image series are converted to through the steadystate spoiled gradient echo (SGPR) signal equation [3], given by
(2) where , , is the repetition time, is the flip angle, is the contrast agent relaxivity taken as 4.2 , is the baseline (precontrast) image intensity, and and are respectively the relaxation and equilibrium longitudinal magnetization that are calculated from a precontrast mapping acquisition.

The undersampled raw (k,t)space data can be related to for a singlecoil data by an undersampling fast Fourier transform (FFT), ,
(3) where represents kspace coordinates.
By simply integrating the three computation steps in (13), we can form a single function modeling the signal evolution in (kt) space given the PK maps , as , where denotes all the predetermined acquisition parameters as mentioned above.
Given the undersampled (k,t)space data , the corrupted image series can be obtained by applying IFFT to , i.e. . We further define a new function that integrates only the first two computation steps (12) to compute the dynamic DCE image series. We will incorporate in our custom loss function that will be explained in the following section.
2.3 PK Parameter Inference via Forward Physical Model Loss
2.3.1 Formulation.
We hypothesize that a direct inversion between corrupted PK parameter maps and is available through forward model, i.e., . However, this cannot provide yet sufficiently accurate estimate of target parameter maps obtained from fullysampled data . To this end, we estimate a correction or residual map from the available signal satisfying . As shown in Fig. 2(a), we observe that residual PK maps involve more sparse representations and exhibit spatially less varying structures inside the brain. The task of learning a residual mapping was shown to be much easier and effective than the original mapping [10]. Following the same approach, we adopt the residual learning strategy using deep CNNs. Our CNN is trained to learn a mapping between and to output an estimate of residual maps ; , where represents the forward mapping of the CNN parameterised by . The final parameter estimate is obtained via .
2.3.2 Loss Function.
We simultaneously seek the signal belonging to the corrected model estimates to be sufficiently close to true signal, i.e., . Therefore, we design a custom loss function which requires solving the forward model in every iteration of the network training. We refer the resulting loss as forward physical model loss. Given a set of training samples of inputoutput pairs (), we train a CNN model that minimizes the following loss,
(4) 
where is a regularization parameter balancing the tradeoff between the fidelity of the parameter and signal reconstruction. We emphasize that the second term in (4) allows the network to intrinsically exploit the underlying contrast agent kinetics in training phase.
2.3.3 Network Architecture.
Figure 3 illustrates our network architecture. The network takes a imagetime series as input, where time frames are stacked as input channels. The first convolutional layer applies filters to each channel individually to extract lowlevel temporal features which are aggregated over frames via learned filter weights to produce a single output per voxel. Following the first layer, inspired by the work on brain segmentation [5], our network consists of parallel dual pathways to efficiently capture multiscale information. The local pathway at the top focuses on extracting details from the local vicinity while the global pathway at the bottom is designed to incorporate more contextual global information. The global pathway consists of dilated convolutional layers with dilation factors of , implying increased receptive field sizes. The filter size of each convolutional layer including dilated convolutions is , and the rectified linear units (ReLU) activation is applied after each convolution. Local and global pathways are then concatenated to form a multiscale feature set. Following this, fullyconnected layers are used to determine the best possible feature combination that can accurately map the input to output of the network. Finally, the last layer outputs the estimated residual maps.
3 Experiments and Results
3.0.1 Datasets.
We perform experiments on fullysampled DCEMRI datasets acquired from three mild ischaemic stroke patients. DCE image series were acquired using a 1.5T clinical scanner with a 3D T1W spoiled gradient echo sequence (TR/TE = ms, flip angle = , FOV = cm, matrix = , slice thickness = mm, sec temporal resolution, dynamics). An intravenous bolus injection of mmol/kg of gadoterate meglumine (GdDOTA) was administered simultaneously. The total acquisition time for DCEMRI was approximately minutes. Two precontrast acquisitions were carried out at flip angles of and to calculate precontrast longitudinal relaxation times.
3.0.2 Preprocessing.
Undersampling was retrospectively applied to the fullysampled data in the  plane using a randomized goldenangle sampling pattern [12] over time (see Fig. 2(b)) with a 10fold undersampling factor. The precontrast first frame was fully sampled. Due to the low temporal resolution of our data, we estimated subjectspecific vascular input functions (VIFs) extracted by averaging a few voxels located on the superior sagittal sinus where the inflow artefact was reduced compared to a feeding artery [4]. Data augmentation was employed by applying rigid transformations on image slices. We generated random 2D+t undersampling masks to be applied on the images of different orientations. This allows the network to learn diverse patterns of aliasing artifacts. All the subject’s data required for network training/testing were divided into nonoverlapping 3D blocks of size , resulting in 64 blocks per subject.
3.0.3 Experimental setup.
All experiments were performed in a leaveonesubjectout fashion. The networks were trained using the Adam optimizer with a learning rate of (using a decay rate of ) for 300 epochs and minibatch size of 4. To demonstrate the advantage of the proposed method, we compare it with the stateoftheart modelbased iterative parameter reconstruction method using the MATLAB implementation provided by the authors [3]. We use the concordance correlation coefficient (CCC) and structured similarity metric (SSIM) metrics to quantitatively assess the PK parameter reconstruction, and peak signaltonoise ratio (PSNR) metric to assess the image reconstruction. Experiments were run on a NVIDIA GeForce Titan Xp GPU with 12 GB RAM.
3.0.4 Results.
Figure 4 shows the qualitative PK parameter reconstructions obtained from different methods using 10fold undersampling. The results indicate that CNN incorporating two loss terms simultaneously produces better maps and considerably higher SSIM score calculated with respect to fullysampled PK maps. The modelbased iterative reconstruction yields the PK maps where the artifacts caused by undersampling are still observable.
In Fig. 5 we present the exemplary reconstructed images obtained by applying the operation to the estimated PK maps. All the reconstruction approaches result in high quality images, however, the modelbased reconstruction can better preserve the finer details. Unfortunately, our fullysampled data suffer from Gibbs artifacts appearing as multiple parallel lines throughout the image. As marked by white arrows, our CNN method can significantly suppress these artifacts whereas they still appear in the image obtained by modelbased iterative reconstruction. Finally, Fig. 8 demonstrates the quantitative results of parameter estimation and image reconstruction. The highest CCC and SSIM values for parameter estimation are achieved by our CNN model when both loss terms are incorporated with and , yielding an average score of 0.88 and 0.92, respectively. The difference is statistically significant for both CCC () and SSIM () when compared against modelbased reconstruction. The modelbased reconstruction performs the highest PSNR for image reconstruction, where it is followed by the proposed CNN with . The difference between them is statistically significant with . The PSNR also shows a decreasing trend with increasing as expected.
We emphasize that the parameter inference of our method on a 3D test volume takes around 1.5 seconds while the modelbased method requires around 95 minutes to reconstruct the same volume, enabling faster computation.
4 Conclusion
We present a novel deep learning based framework for direct estimation of PK parameter maps from undersampled DCE imagetime series. Specifically, we design a forward physical model loss function through which we exploit the physical model relating the contrast agent kinetics to the timeresolved DCE signals. Moreover, we utilize the residual learning strategy in our problem formulation. The experiments demonstrate that our proposed method can outperform the stateoftheart modelbased reconstruction method, and allow almost instantaneous inference of the PK parameters in the clinical workflow of DCEMRI.
4.0.1 Acknowledgements.
The research leading to these results has received funding from the European Unions H2020 Framework Programme (H2020MSCAITN 2014) under grant agreement no 642685 MacSeNet. We also acknowledge Wellcome Trust (Grant 088134/Z/09/A) for recruitment and MRI scanning costs.
References
 [1] Fang, R., et al.: Direct estimation of permeability maps for lowdose ct perfusion. In: IEEE ISBI. pp. 739–742 (April 2016)
 [2] Guo, Y., et al.: Highresolution wholebrain DCEMRI using constrained reconstruction. Med. Phys. 43(5), 2013–2023 (2016)
 [3] Guo, Y., et al.: Direct estimation of tracerkinetic parameter maps from highly undersampled brain dynamic contrast enhanced MRI. MRM 78(4) (2017)
 [4] Heye, A.K., et al.: Tracer kinetic modelling for dcemri quantification of subtle bloodâbrain barrier permeability. NeuroImage 125, 446 – 455 (2016)
 [5] Kamnitsas, K., et al.: Efficient multiscale 3d CNN with fully connected CRF for accurate brain lesion segmentation. Medical Image Analysis 36, 61–78 (2017)
 [6] Lebel, R.M., et al.: Highly accelerated dynamic contrast enhanced imaging. MRM 71(2), 635–644 (2014)
 [7] O’Connor, J.P.B., et al.: Dynamic contrastenhanced MRI in clinical trials of antivascular therapies. Nat. Rev. Clin. Oncol. 9(3), 167–77 (Mar 2012)
 [8] Patlak, C.S., et al.: Graphical evaluation of bloodtobrain transfer constants from multipletime uptake data. J. Cereb. Blood Flow Metab. 3(1), 1–7 (1983)
 [9] Sourbron, S.P., Buckley, D.L.: Classic models for dynamic contrastenhanced MRI. NMR in Biomedicine 26(8), 1004–1027 (2013)
 [10] Zhang, K., et al.: Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising. IEEE Trans. Image Process. 26(7), 3142–3155 (July 2017)
 [11] Zhang, T., et al.: Fast pediatric 3d freebreathing abdominal dynamic contrast enhanced mri with high spatiotemporal resolution. JMRI 41(2), 460–473 (2015)
 [12] Zhu, Y., et al.: GOCART: GOldenangle CArtesian randomized timeresolved 3D MRI. Magn. Reson. Imag. 34(7), 940 – 950 (2016)