Epicardial adipose tissue (EAT), the fat located between the myocardium and the visceral pericardium [1] which serves as an energy store [2], has been hypothesized as a contributor to the inflammatory burden via paracrine mechanisms [3]. EAT has been suggested as an independent marker for cardiovascular risk [2, 4, 5, 6, 7, 8] and, in particular, as an independent atrial fibrillation (AF) prognostic marker [9]. In computed tomography (CT), EAT volume (EATv) and density (EATd) are parameters that are often used to quantify EAT [8]. EATv refers to the extent of EAT accumulation and increased EATv has been found to correlate with the prevalence and the recurrence of AF after ablation therapy [1, 10, 11, 12, 13, 14, 15]. In addition to AF, increased EATv is also associated with atherosclerosis [16, 17], carotid stiffness [18], myocardial infarction [19], and coronary artery calcification [20, 21]. EATv has been associated with severity of coronary artery disease (CAD) [6, 22, 23]. A higher EATd (i.e., radiodensity of EAT) in CT images is correlated with inflammation as a result of arrest of lipid maturation [24]. Researchers have also suggested a link between EATd and high risk mortality and plaque presence and progression [24].

These health risk factors emphasize the need for direct quantification of EAT (i.e., EATv and EATd). However, manual quantification is time-consuming to accomplish in clinical practice in the light of the high workload on physicians and radiographers, and so EAT is not routinely quantified. Automation of the quantification task diminishes the variability in readings introduced by different observers and removes the high dependence of state-of-the-art methods on user interaction for EAT segmentation, resulting in high reproducibility of studies and less time-consuming analysis. In general, fully automated quantification of EATv and EATd requires advanced techniques and, in this work, we propose a deep learning (DL) framework for carrying out the estimation of these two quantities autonomously. DL techniques, a set of machine learning methods, have proved to be very effective for automated detection and segmentation of a wide range of medical images with a high degree of accuracy [25, 26, 27, 28, 29]. The proposed methodology, therefore, mainly consists of image classification and segmentation DL-based models to perform the desired quantification. We compare the performance of our approach to other DL approaches proposed in the literature for accomplishing EAT quantification.

The segmentation of EAT in cardiac CT (CCT) image slices is important for EAT quantification and various semi-automatic segmentation approaches have been developed [30, 31, 32]. An overview of these approaches can be described as follows: after initial preprocessing involving the removal of all other structures in the CT images apart from the heart, these methods require an expert to scroll through the CT slices to identify some control points along the border of the pericardium, then use some interpolation methods (such as the cubic spline function techniques) to obtain smooth pericardial contour, and then identify the pericardial fat by thresholding.

Fully automated and semi-automated non-DL based EAT segmentation approaches have also been developed [33, 34, 35, 36, 37, 38]. The time taken for obtaining the segmentation masks per patient, according to [35], could be more than 15 minutes for such fully automated approaches. The DL-based automatic EAT quantification methods that have been reported for EAT segmentation or quantification include those by Commandeur et al. [39] and Li et al. [40]. In Commandeur et al. [39], two DL-based models were developed; one is used to determine heart limits and perform heart segmentation (i.e., thoracic mask segmentation) and the other was used in combination with a statistical shape model for the detection of the pericardium. EAT was then quantified by further post-processing using thresholding [–190, –30] HU. In Li et al. [40], a DL-based model was developed for the segmentation of the pericardium across multiple adjacent slices using multiple slices as input to the model. A smoothing operation is then employed by finding a solution to a partial differential equation of the 3-dimensional gradient vector flow in order to reduce the prediction of false positive and negative regions in the segmented pericardial images. EAT is then deduced by thresholding [–175, –15] HU.

To our knowledge, there is no unified method in the literature capable of both autonomous EATv and EATd. Methods exist for the EAT segmentation and EATv estimation [39, 40, 41, 42], but our approach differs from these methods in that EAT segmentation does not require any further post-processing (e.g., thresholding or smoothening operations with filters) after the prediction with the DL-based EAT segmentation model. Also, no previous work has attempted to estimate EATd alone with DL nor combine the estimation of EATv and EATd. Although some authors, such as [41, 42], have stated their quantification of EATd as the mean attenuation of EAT segmented using DL models, we are not aware of any previous work that demonstrates or carried out practical analysis of the EATd quantification using DL. Our results in the analysis of EATd strengthens the correctness of our approach of EAT segmentation and emphasizes the correctness of the results obtained for EATv estimation.

We consider the ResNet (ResNet50) [48] DL architecture for the image classification task of selecting images containing the EAT from the set of all CT images acquired for a given patient, eliminating slices above the bifurcation of the pulmonary trunk or below the cardiac apex. That is, the model classification automatically identifies slices above the superior extent of the left main coronary artery and also those slices below the cardiac apex for elimination in further analysis. A detailed description of the ResNet50 architecture is given in Supplementary Table 1 and Supplementary Materials.

During model training, the input CCT images were resized to 224 $\times$ 224 and, as part of the in-training data augmentation, we rotated the images up to $\pm$15° and their intensities normalized. The weights of the models were randomly initialized. Training was carried out for 60 complete epochs using a batch size of 128 images. The binary cross-entropy function was used as the loss function. The Root Mean Squared Propagation (RMSProp) optimizer was chosen as the optimisation technique and set to an initial learning rate of 1e-4 with a step decay schedule.

For our experiments in this study, our 300-patient dataset is divided into two subsets, Dataset 1 (41,979 CT slices) and Dataset 2 (36,428 CT slices), with each subset consisting of 150 patients. We used Dataset 1 for training, validating and evaluating the DL models, and Dataset 2 for evaluating EATv and EATd estimations. Splitting our dataset in this way ensures that in the estimation of EATv and EATd, the DL models trained on Dataset 1 have not seen the CT slices in Dataset 2 during training, allowing us to have a proper evaluation of the proposed methods for EATv and EATd quantification.

The CT image dataset of 150 patients (Dataset 1) of which 23,771 images contain the EAT were used to train, validate and test the ResNet50 model. In particular, 15% of the images were selected; these were then divided into two equal sets representing the validation and test sets. Thus, the number of the training, validation and test (evaluation) images are 35,683, 3148 and 3148, respectively.

Of the 35,683 images in the training set, the number of images with the EAT present and absent are 20,234 and 15449, respectively. We used weighted loss function to address this data imbalance. Let $\{(x_1,y_1),(x_2,y_2),\mathrm{\dots },(x_n,y_n),\mathrm{\dots },(x_N,y_N)\}$, where $x$ is the two-dimensional input image, denote a training set of $N$ samples and $y\in {\{0,1\}}^C$ is a binary one-hot encoded label (i.e., $C$ = 2 in the present case), then the weighted loss function can be written as follows:

(1)$$\begin{array}{c}\hfill E_w\left(\theta \right)=-\frac{1}{N}[\lambda {}_0\sum {}_{n=1}{}^{N}T_0(x_n)y_{n}log({\widehat{y}}_n(x_n,\theta ))\hfill \\ \hfill +{\lambda }_1\sum {}_{n=1}{}^{N}T_1(x_n)y_{n}log({\widehat{y}}_n(x_n,\theta ))]\hfill \end{array}$$

where $\theta$ denotes the trainable model weights; ${\widehat{y}}_n(x_n,\theta )$ is the posterior probability obtained after applying the sigmoid activation function on the model’s output layer; $T_0(x_n)$ and $T_1(x_n)$ are functions indicating whether image $x_i$ belongs to class 0 or class 1, respectively (i.e., whether the EAT is absent or present in image $x_i$, respectively); and ${\lambda }_0$ and ${\lambda }_1$ are weights for penalizing the loss function for false negatives and false positive errors, respectively (i.e., $x_i$ that belongs to class 0 is wrongly classified as belonging to class 1 or belongs to class 1 is wrongly classified as belonging to class 0, respectively). The weights, ${\lambda }_0$and ${\lambda }_1$, are given as follows:

(2)$${\lambda }_i=\frac{1}{k_i}\cdot \frac{N}{C}$$

where $k_i$ denotes the number of images in class $i$. In our case, $N=35952$; class 0 and class 1 are subgroups indicating the collection of samples where EAT is absent and present, respectively; then, ${\lambda }_0=(1/20234)\times (35683/2)=0.882$ and ${\lambda }_1=(1/15449)\times (35683/2)=1.155$. That is, the images contain the EAT (class 1) are weighted as being more significant than those without the EAT (class 0).

The image segmentation task of masking the regions of the EAT from the CT slices is performed by the UNet DL model [49]. A detailed description of the UNet architecture is given in Supplementary Table 2 and Supplementary Materials. Briefly, the architecture includes batch normalization to enhance robustness of the model [50] and ‘dropout’ operation [51] to avoid problems associated with model overfitting.

During model training, the input CCT images were resized to 224 $\times$ 224 and image rotation, intensity normalization and cropping were performed as part of the in-training data augmentation. The model parameters were randomly initialized, and training proceeded for 60 epochs using a batch size of 64 images. The sparse categorical cross-entropy function and the Adam optimizer were used as the loss function and the optimisation method, respectively, with an initial learning rate of 1e-4 decreasing exponentially at a rate of –0.05 after the first 5 epochs.

For $N$ samples of the training set $\{(x_1,y_1),(x_2,y_2),\mathrm{\dots },(x_n,y_n),\mathrm{\dots },(x_N,y_N)\}$ where $x$ is the two-dimensional input image and $y$ is the two-dimensional segmentation mask, the sparse categorical cross-entropy loss function can be written as follows:

(3)$$E_w\left(\theta \right)=-\frac{1}{N}\left[{\sum }_{n=1}^N{y}_n\log ({\widehat{y}}_n(x_n,\theta ))\right]$$

where $\theta$ denotes the trainable parameters of the model and ${\widehat{y}}_n(x_n,\theta )$ is the posterior probability obtained following ‘sigmoid’ activation function on the output layer of the model.

The CT image dataset of 150 patients in Dataset 1 with the EAT present (23,771) were used to train, validate and test the EAT segmentation models. In particular, 15% of these images were selected and divided into two equal sets namely, the validation and test sets. Thus, the numbers of images in training, validation and test sets are 20,207, 1782 and 1782, respectively.

The dice score, or dice similarity coefficient (DSC), is a measure of similarity between the label and predicted segmentation masks. DSC score can be written as follows:

(4)$$DSC=\frac{2\left|A\cap B\right|}{\left(\left|A\right|+\left|B\right|\right)}$$

where $A$ and $B$ represent the two sets (images), $\left|A\right|$ and $\left|B\right|$ represent the cardinalities (i.e., the number of elements) of set $A$ and $B$, respectively. The DSC is a useful measure of spatial overlap often used in image segmentation to quantify the accuracy of the predicted mask with respect to the ground truth mask. It is also used as a statistical validation metric by computing the DSC of several images obtained using images from evaluation subset and computing the mean DSC.

The estimation of EATv involves the integration (summation) of the interslice volumes. Each interslice volume is approximated by computing the volume between two consecutive slices using the following equation:

(5)$$V={\sum }_i^N{v}_i$$

where $V$ is the estimated EATv and $N$ is the number of slices. $v_i$ (the $i$th interslice volume, i.e., $v_i$ the arithmetic mean of the areas of two consecutive slices multiplied by the distance between them in the direction of the $z$-axis) is computed as follows:

(6)$$v_i=\frac{A_i+A_{i+1}}{2}{z}_{k_i}$$

$z_{k_i}$ is the distance between two consecutive slices in the $z$ direction and can be expressed as $z_{k_i}=z_i\widehat{k}$, where $\widehat{k}$ is the unit vector in the $z$ direction. Also, $z_i$ is the perpendicular distance between two consecutive parallel slices; $A_i$ and $A_{i+1}$ represent the areas of the two consecutive slices, and $A_i$ is defined as follows:

(7)$$A_i=n_i{s}_x{s}_y$$

where $n_i$ is the number of pixels that constitutes the EAT area on the $i$th slice; $s_x$ and $s_y$ are the pixel dimensions in the $x$ and $y$ directions, respectively. If $s=s_x=s_y$, then $A_i=n_i{s}^2$.

EATd is estimated by computing the mean attenuation of EAT across all the slices. The intensity values of the EAT for each slice, with range [–190, –30], is extracted using its segmentation mask obtained from the EAT segmentation model to mark the ROI. EATd can be expressed as follows:

(8)$$EATd=\frac{1}{M}{\sum }_{i=1}^{Ns}{\sum }_{j=1}^{N{p}_i}{x}_i(j)$$

where $x_i(j)$ is the intensity value of pixel $j$ in the ROI of slice $i$; $N{p}_i$ is the number of pixels in the ROI of slice $i$; $Ns$ is the number of slices; and $M$ is the total number of pixels for the of the totality of EATv and can be expressed as:

(9)$$M={\sum }_{i=1}^{Ns}N{p}_i$$

Since the size of the input image and output segmentation mask of the segmentation model is 224 $\times$ 224, the output segmentation mask was resized to the size of the CT image (512 $\times$ 512) using area interpolation (which, in this case of image enlargement, is a bilinear interpolation image processing technique involving resampling using pixel area relation [52]) before extracting the intensity values of the pixels within the ROI on the CT image.

The performance metrics of the ResNet50 classification model on the evaluation dataset of N = 3148 are precision (0.980), recall (0.986) and F${}_1$ score (0.983) and the confusion matrix is given in Table 2. Some examples of the predictions of the classification model are given in Fig. 2.

Fig. 2.

Some examples of the prediction of the classification model. Two examples are given for each of the cases (0 – absence of EAT, 1 – presence of EAT): (0/0) (i.e., ground truth/prediction), (0/1), (1/0) and (1/1). For the (1/0) and (1/1) cases, the images to the right show the EAT in red as given by an expert human reader. The (0/1) and (1/0) cases are images which the classification model got wrong.

The performance metrics of the segmentation model on the 1782 evaluation set are given as follows: the mean DSC is 0.844 ($\pm$ 0.19 standard deviation); 25%, 50% (median) and 75% percentiles of DSC are 0.81, 0.84 and 0.87, respectively; the maximum DSC is 0.95. To ensure that the proposed segmentation model is robust, the mean ($\pm$ std.), maximum and median DSC of this model (0.844 ($\pm$ 0.19), 0.95 and 0.84, respectively) is compared with the results from a 5-fold cross-validation of the same model architecture are presented in Table 3 (i.e., 20% of Dataset 1 is used as subsamples/validation set and each of the 5 subsamples used exactly once, creating 5 models); the best of the five models (model 4) gives the mean ($\pm$ std.), maximum and median DSC as 0.851 ($\pm$ 0.17), 0.96, and 0.84, respectively. Fig. 3 gives some examples of the predictions of the segmentation model with varying DSC.

Fig. 3.

Some examples of the predictions of the segmentation model. The corresponding dice scores are shown at the bottom of each of the examples.

The regression, the kernel density estimates and the Bland-Altman plots of the predicted volumes against the label volumes computed for the study population of Dataset 1 (the ‘training’ dataset) and Dataset 2 (the ‘evaluation’ dataset), each of which consists of 150 patients, using the proposed framework are shown in Fig. 4. Fig. 4a,c show the regression plots with volume estimated using Eqns. 5,6,7. In our case, $s_x=s_y$ and $z_{k_i}$ is approximately between 2.2 mm and 4.5 mm. The Pearson correlation coefficient, $\rho$, between the label and predicted EATv are 0.977 (R${}^2$ = 0.954) and 0.971 (R${}^2$ = 0.943) for the populations in Dataset 1 and Dataset 2, respectively. The arithmetic means (93.36 mL for label and 101.16 mL for prediction in Dataset 1, and 105.30 mL for label and 103.29 mL for prediction in Dataset 2) of the distributions (the dashed vertical lines on the kernel density estimates plot of Fig. 4b,e) show close agreements between the label and prediction of the mean EATv of the populations. Bland-Altman plots Fig. 4c,f indicates 95% of the difference between the label and predicted EATv values will be between –22.67 mL and +17.06 mL of the mean EATv difference value of –2.81 mL for Dataset 1, and between –21.13 mL and +25.15 mL of the mean EATv difference value of +2.01 mL for Dataset 2.

Fig. 4.

The plots of the predicted volume against the label volume. Plot (a) represents the regression plot; plot (b) represents the kernel density estimates and histogram plots of the two variables (label volume and predicted volume) with the dashed vertical lines representing the arithmetic mean of the distributions. The symbols $\rho$ and $p$ represent the p-value and Pearson correlation coefficient, respectively. Plot (c) represents the Bland-Altman plot of the predicted volumes against the label volumes where the lower and upper dashed horizontal lines are the confidence interval at 95%. (a–c) show the plots for Dataset 1. For Dataset 2, plot (d) represents the regression plot; plot (e) represents the kernel density estimates and histogram plots of the two variables (label volume and predicted volume); plot (f) represents the Bland-Altman plot of the predicted volumes against the label volumes.

Similarly, the regression, the kernel density estimates and the Bland-Altman plots of the predicted EATd versus the label EATd estimated for the study population of Dataset 1 and Dataset 2 are given in Fig. 5. Fig. 5a,c show the regression plots, and the Pearson correlation coefficient between the label and predicted EAT mean attenuation are 0.964 (R${}^2$ = 0.930) and 0.972 (R${}^2$ = 0.945) for the populations in Dataset 1 and Dataset 2, respectively. The arithmetic means (–84.51 HU for label and –85.07 HU for prediction in Dataset 1, and –85.78 HU for label and –86.44 HU for prediction in Dataset 2) of the distributions (in Fig. 5b,e) show close agreements between the label and prediction of the mean attenuation of the populations. Bland-Altman plots Fig. 5c,f) indicates 95% of the difference between the label and predicted EATd values will be between –3.13 HU and +4.25 HU of the mean EATd difference value of 0.56 HU for Dataset 1, and between –2.31 HU and +3.36 HU of the mean EATd difference value of 0.66 HU for Dataset 2. Moreover, the estimation of EATv and EATd took approximately 39.54 sec and 30.74 sec per patient, respectively. Some examples of the label (ground truth) EATv and EATd values and the predicted values from the proposed framework are given in Supplementary Table 3.

Fig. 5.

The plots of the predicted attenuation against the label attenuation. Plot (a) represents the regression plot and plot (b) represents the kernel density estimates and histogram plots of the two variables (label and predicted mean attenuations) with the dashed vertical lines representing the arithmetic mean of the distributions. The symbols $\rho$ and $p$ represent the p-value and Pearson correlation coefficient, respectively. Plot (c) represents the Bland-Altman plot of the predicted versus the label mean attenuations where the lower and upper dashed horizontal lines are the confidence interval at 95%. (a–c) show the plots for Dataset 1. For Dataset 2, plot (d) represents the regression plot; plot (e) represents the kernel density estimates and histogram plots of the two variables (label volume and predicted volume); plot (f) represents the Bland-Altman plot of the predicted volumes against the label volumes.

The performance of the method proposed in [39] was evaluated for EAT segmentation on 10% of the dataset of 250 subjects (i.e., 25 in the evaluation set) and gave a median DSC of 0.823. EATv on the 250 patients gave a correlation of (R${}^2$0.924). The 95% limits of agreement ranged from Bland-Altman plots is approximately –27 mL to 23 mL. The approach presented in this paper significantly differs with the method proposed in [39] in that, unlike in [39] where the task of quantifying EAT is divided into 3 separate tasks (namely, slice selection task, heart localization task, and pericardium line detection) and a further post-processing of the outcome via thresholding, our approach only involved two tasks: slice selection and EAT segmentation. To emphasize, our approach does not require any further post-processing. We note that the goal of the slice selection task in our method and the method in [39] are the same except that we have used a more recent classification architecture (i.e., ResNet50 [48]). Moreover, the authors of [39] have extended their approach to a larger multi-centre cohort in [41] and have reported a median DSC of 0.873 for EAT segmentation on 10% of 614 CT dataset and EATv evaluation of (R${}^2$0.974) on a dataset of 614 studies. The 95% limits of agreement ranged from Bland-Altman plots is –19.59 mL to 21.42 mL. A further extension of the work in [39] is given [42].

Also, the method proposed in [40] was trained on a dataset of 88 subjects, and reported a mean DSC score of 0.973 using an evaluation dataset of only 15 subjects; thus, the evaluation dataset is not large enough for making a fair comparison with other methods. In addition, our approach differs from the method presented in [40] in that it does not require smoothing operation by solving a differential equation nor any post-processing step via thresholding.

Our models were trained on a dataset from a single hospital system. The data augmentation techniques, the batch normalization and dropout operations and the weighted loss function for addressing data imbalance used during model training may have improved the chance that performance may not deteriorate significantly from datasets from elsewhere but it would be useful assessing the performance of the models using an external dataset. If needed, the performance of the models may then be improved with training on multi-centre datasets to enhance generalizability with little or no modification to the proposed methods. Approaches that may be explored to address model generalization issues include transfer learning and federated-learning [53].

In relation to gender, men constitute the majority (66%) of the dataset we have used for model training. In the evaluation of our models trained with this dataset using 5-fold cross-validation (Table 3), there is no indication of biasness of the models at slice level towards a particular gender. Biases in the training datasets can affect performance of DL models. As such, future work on this research would focus on investigating the estimation of EATv and EATd at patient-level for biasness. A possible approach for addressing biasness in CT images is by balancing the dataset using generative models in the form of data augmentation. An in-depth discussion on generative models is beyond the scope of this paper and we refer readers to [54] for more details on this technique.

Our framework focused on EATv and EATd quantification for CT images, future work will focus on training using heterogenous multi-centre dataset as well as on analysis of EAT for cardiovascular risk and outcome prediction. Future direction of this work will also include using machine learning methods for quantifying the distribution of EAT given that the location of EAT is a disease-specific risk factor (e.g., thickness of peri-atrial EAT being a predictor of AF recurrence [55]).