3.1.

Patient Population

The cohort used in this study comprised 264 patients, who underwent chest CT and had a positive RT-PCR test result for SARS-CoV-2. A total of 197 patients were included in the training cohort ($N_{\mathrm{cov}}$), and 68 were used for external validation ($N_{\mathrm{ext}}$). Datasets for 187 out of 197 patients from the training cohort were collected from the prospective, international, multicenter registry involving centers from North America [Cedars-Sinai Medical Center, Los Angeles ($n=75$)], Europe [Centro Cardiologico Monzino ($n=64$), and Istituto Auxologico Italiano ($n=17$); both Milan, Italy], Australia [Monash Medical Centre, Victoria, Australia ($n=6$)], and Asia [Showa Medical University, Tokyo, Japan ($n=25$)], where either non-contrast ($n=157$) or contrast-enhanced ($n=30$) chest CT was performed to aid in the triage of patients with a high clinical suspicion for COVID-19, in the setting of a pending RTPCR test or comorbidities associated with severe illness from COVID-19. The population is given in Table 1. Datasets for the remaining 10 COVID-19 patients were derived from an open-access repository of non-contrast CT images; therefore, no clinical data were provided for this cohort. Out of 31,560 transverse slices available, 15,588 had lesions. The external testing cohort comprised 68 non-contrast CT scans of COVID-19 patients: about 50 from an open-access repository²⁹ and 18 additional ones from Italy (Centro Cardiologico Monzino). There were 12,102 transverse slices available in this cohort, and 6,503 had lesions (Table 2). All data were deidentified prior to being enrolled in this study. The CT images from each patient and the clinical database were fully anonymized and transferred to one coordinating center for core lab analysis. The study was conducted with the approval of local institutional review boards (Cedars-Sinai Medical Center IRB# study 617) and written informed consent was waived for fully anonymized data analysis.

Table 1

Patient baseline characteristics and imaging data in a training cohort.

Table 2

Image findings.

3.2.

Ground Truth Generation

Images were analyzed at the Cedars-Sinai Medical Center core laboratory by two physicians (K.G. and A.L.) with 3 and 8 years of experience in chest CT, respectively, and who were blinded to clinical data. A standard lung window (width of 1500 HU and level of $-400\mathrm{HU}$) was used. Lung abnormalities were segmented using semi-automated research software (FusionQuant Lung v1.0, Cedars-Sinai Medical Center, Los Angeles, California). These included GGO, consolidation, or pleural effusion according to the Fleischner Society lexicon. Consolidation and pleural effusion were collectively segmented as high-opacity to facilitate the training of the model due to a limited number of slices involving these lesions. Chronic lung abnormalities, such as emphysema or fibrosis, were excluded from segmentation, based on correlation with previous imaging and/or a consensus reading. GGO was defined as hazy opacities that did not obscure the underlying bronchial structures or pulmonary vessels; consolidation as opacification obscuring the underlying bronchial structures or pulmonary vessels; and pleural effusion as a fluid collection in the pleural cavity. The total pneumonia volume was calculated by summing the volumes of the GGO and consolidation components. The total pneumonia burden was calculated as (total pneumonia volume/total lung volume) × 100%. Difficult cases of quantitative analysis were resolved by consensus.

3.3.

Controls Dataset

Additionally, to assess the diagnostic performance of the methods trained and tested with controls (without any lung abnormalities), we utilized a set of $N_{\text{control}}=695$ cases from the national lung screen trial (NLST)³⁰ with normal lung scans. The population characteristics are described in Table 3.

Table 3

NLST controls baseline characteristics.

4.1.

Data Preprocessing

CT scans from different scanners or with different reconstruction parameters may have different appearance (as seen in column 1 of Fig. 1) and contain voxel values (HU) ranging between $-1024$ to $+3071$ for a 12-bit scan. Therefore, there is a need for homogenizing the data before we train or infer from it. The input stack of CT images $\mathbf{I}$ are first cropped to the body region of the middle-most image and resized to $256\times 256$. Because we have a very small dataset to train on, we randomly augment the data with rotation of $[-10\mathrm{deg},+10\mathrm{deg}]$, translation of up to 10-pixels in the $x$- and $y$-directions, and scaling of [0.9, 1.05] times. Finally, we normalize the data by clipping the Hounsfield units between $-1024$ to $+600$ (expert reader’s lung window), followed by a voxel intensity scaling technique called standardization or $Z$-score normalization.

Fig. 1

Order of data preprocessing from input $\mathbf{I}$ (left) to processed output ${\mathbf{I}}_{\mathrm{std}}$ (right).

To crop the scan to the body region, we threshold the scan at $-500\mathrm{HU}$ and create a binary mask followed by a series of morphological operations: closing, erosion, dilation, etc., and obtain a bounding box around the largest object in the threshold scan. Transferring this bounding box, the original scan returns the body cropped input scan (shown in column 2 of Fig. 1).

4.2.

Network Architecture

The network architecture, shown in Fig. 2, is inspired by the hierarchical multi-scale attention for semantic segmentation²⁶ with major changes in the attention branch. Instead of the attention branch looking at the input at various different scales as in Ref. 26, we formulate the attention branch to focus on adjacent slices to aggregate information about the lesions/anatomy from the neighboring slices using a ConvLSTM in the attention branch of the network to improve lesion recognition.

Fig. 2

Framework of the proposed method.

4.2.2.

Attention branch ${\mathbf{\Phi }}_{\mathrm{attn}}$

All of the errors made by the main branch in ambiguous and difficult to segment parts of the lesions are corrected by the attention branch using information from adjacent slices (shown in Fig. 3). The attention branch comprises a sequential processor ${\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{clstm}}$, a segmentation block ${\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{seg}}$, and a self-attention block ${\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{attn}}$.

Eq. (5)

$$\alpha ={\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{attn}}({\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{clstm}}(\mathbf{I})),$$

where $\alpha$ is the self-attention.

Eq. (6)

$$\begin{gathered} {\mathbf{S}}_{\mathrm{attn}}={\mathrm{\Phi }}_{\mathrm{attn}}(\mathbf{I}), \\ ={\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{seg}}({\mathbf{I}}_{\mathbf{k}}\oslash {\mathrm{\Phi }}_{\mathrm{attn}}^{\mathrm{clstm}}(\mathbf{I})), \end{gathered}$$

where ${\mathbf{S}}_{\mathrm{attn}}\in {\mathbb{R}}^{H\times W\times C}$ are the output features from the segmentation block 2 and $C$ is the number output classes.

Fig. 3

Intermediate output showing error correction by the attention branch for four different cases in each row. Blue indicates GGOs, and yellow indicates high opacities. Errors are encircled in red.

Convolutional LSTM ${\mathbf{\Phi }}_{\mathrm{attn}}^{\mathrm{clstm}}$

We used ConvLSTM³³ for processing sequential data. The ConvLSTM block allows for imitating a radiologist reviewing adjacent slices of a CT scan and aggregate lesion information from adjacent slices to detect lung abnormalities and ensure appropriate annotations.

Segmentation block 2 ${\mathbf{\Phi }}_{\mathrm{attn}}^{\mathrm{seg}}$

This block is structurally identical to segmentation block 1, except for the input layer. It takes in the main segmentation slice concatenated with ConvLSTM output as the input.

Attention block ${\mathbf{\Phi }}_{\mathrm{attn}}^{\mathrm{attn}}$

As in Ref. 26, we also adopt an attention mechanism to combine multi-branch outputs (${\mathbf{S}}_{\text{main}}$ and ${\mathbf{S}}_{\mathrm{attn}}$) together at a pixel level. The attention block is identical to the segmentation block in structure with the only difference being that the final $1\times 1$ convolutional layer is followed by a sigmoid layer. This block takes in the output of the ConvLSTM block, as shown in Eq. (5), as input and learns to pixel-wise weight ($\alpha$) the outputs from the two branches to produce the final prediction [Eq. (7)].

The final prediction is given by the following equation in which the argmax is taken over the channel dimension:

Eq. (7)

$${\mathbf{S}}_{\text{out}}=\underset{C}{\arg \max }(\sigma ((1-\alpha ){\mathbf{S}}_{\text{main}}+\alpha {\mathbf{S}}_{\mathrm{attn}})),$$

where $\sigma :{\mathbb{R}}^C\to {(\mathrm{0,1})}^C$ is the Softmax over the channel dimension.

4.3.

Loss Function

In our training, we utilize a combination of focal loss³¹ and Visual Geometry Group (VGG) loss.³⁴ The focal loss compensates for the imbalance between background, GGO, and high opacity classes. The importance for each of the classes in focal loss was set to [0.1, 1.0, 1.0], respectively, and the focusing parameter $\gamma$ was set to 3. This focusing parameter in the focal loss allows the model to penalize the hard to classify samples more than the easy ones. We tap into the low-level features in the VGG network to compute the VGG loss, which represent edge information, for better segmentation output. These losses are weighted equally $(\lambda =1.0)$ during training

4.4.

Optimization

The model parameters were optimized using an Adam (adaptive moment estimation) optimizer³⁵ with initial learning rate of ${10}^{-3}$, weight decay of ${10}^{-6}$, and training batch size of 32. All of the model parameters were initialized using Xavier initialization,³⁶ except for the dense block, which was initialized using the weights pre-trained on ImageNet.³⁷ To avoid over-fitting while fully train the model, we use a popular learning rate scheduler called ReduceOnPlateau (Fig. 4).³⁸ In this technique, a metric (validation loss, accuracy, etc.) is continuously monitored throughout the training. If no improvement is seen in the tracked metric for “patience” number of epochs/iterations, the current learning rate is then reduced by the given “factor.” The training continues as usual until the learning rate is reduced beyond a certain minimum (${10}^{-7}$). As soon as the learning rate hits this minimum, the training is stopped, saving the model at the last best validation metric step. In our experiment, the parameters factor and patience were set to 0.1 and 5, respectively.

Fig. 4

Reduce-on-plateau learning rate scheduler. Bad epochs refers to the number of epochs for which the validation loss has not decreased.

4.5.

Implementation

We trained the model using the Pytorch (v1.7.1) deep-learning framework and incorporated research CT lung analysis software (Deep Lung) written in C++. The training was performed on an NVIDIA TITAN RTX 24GB GPU with a tenth generation Intel Core i9 CPU. Deep Lung can be used with or without the GPU acceleration.

5.1.

Five-Fold Cross-Validation

The primary endpoint of this study was the performance of the deep-learning method compared with the evaluation by the expert reader. The model is extensively evaluated using the Dice similarity coefficient for structural similarities. The reported Dice score is the mean of per-patient Dice scores computed over all slices in the scan. We also show the quantitative performance on volumes using the Bland–Altman plot and coefficient of determination $R^2$ (Pearson correlation). To perform a robust non-biased evaluation of the framework, five-fold cross-validation was used, using five independently trained identical models and five exclusive hold-out sets, each of 20%. The whole cohort of $N_{\mathrm{cov}}=197$ cases was split into five subsets called folds. For each fold of the five-fold cross-validation, the following data splits were used: (1) training split (125 or 126 cases) was used to train the ConvLSTM; (2) alidation split (32 cases) was defined to tune the network, select optimal hyperparameters, and verify that there was no over-fitting; and (3) test split (39 or 40 cases) was used for the evaluation of the method. The final results were obtained by concatenating the results from five test subsets. Thus, the overall test population was 197, referred to as internal test set further in the paper. We also test our model on an unseen external dataset consisting of $N_{\mathrm{ext}}=68$ patients.

5.2.

Diagnostic Per-Patient Performance

To assess the diagnostic performance of the convLSTM on a per-patient basis, we trained our model utilizing an additional $N_{\text{control}}^{\text{train}}=197$ NLST controls (read as number of controls in training) during the five-fold cross-validation, making the total training cases $N=N_{\mathrm{cov}}+N_{\text{control}}^{\text{train}}=394$. An additional set of $N_{\text{control}}^{\text{test}}=498$ normal NLST cases (read as number of controls in testing), added during testing, were evaluated with the best fold model from the five-fold cross-validation. Thus, the total normal NLST cases included in experiment sums to $N_{\text{control}}=N_{\text{control}}^{\text{train}}+N_{\text{control}}^{\text{test}}=695$. Each normal case was evaluated with the model, which did not include these cases for training. We report the specificity at 95% sensitivity for the convLSTM models trained with and without additional controls. The diagnostic sensitivity and specificity was compared using McNemar’s test³⁹ on paired measurements.

6.1.

Ablation Study

Table 4 shows how the results are affected by altering different building blocks of our model.⁴⁰ We select the model with the best validation Dice score (mean of GGO and high opacity) for the final evaluation. The model configurations with the highest and lowest performances are highlighted in green and orange, respectively. We experimentally found that the best results were obtained at buffer size $F=3$.

Table 4

Ablation study on fold-1 for model selection (Ncov=197). Highest and lowest performances are highlighted in bold and italic, respectively.

6.2.

Model Comparison

In Table 5, we show the performance of our model as compared with UNet2D and UNet3D across five-folds ($N_{\mathrm{cov}}=197$). For fair comparison, UNet2D and UNet3D were trained with an identical training setup to our model, i.e., the same loss function, optimizer, learning rate strategy, and training fold splits. The performance is measured with two main metrics: mean Dice score and compute resource utilization. The mean Dice score reported in Table 5 gives the binarized mean Dice score per class. The computation time and memory are calculated for 128 CT slices and 16 CT slices, respectively, on an Nvidia Titan RTX GPU and Intel i9 CPU using Pytorch Profiler.⁴¹ In Fig. 5, we show the significance of our results using the Wilcoxon signed-rank test. We see that our model outperforms UNet2D ($p=0.001$) and Unet3d ($p<0.0001$) in segmenting high-opacities, has a comparable performance to UNet2D ($p=0.22$) in segmenting GGOs, and significantly outperforms UNet3D ($p<0.0001$) in segmenting GGOs. But the major advantages of our model over the other two are in terms of computational resources as follows:

Hence, it is can be easily deployed on less powerful machines in clinical setups.

Table 5

Model comparisons on Ncov=197 (UNet2D, UNet3D, and our). Best performance is highlighted in bold.

Fig. 5

Box plot for $N_{\mathrm{cov}}=197$ cases displaying the significance of Dice scores between models using the Wilcoxon signed-rank test for (a) GGO and (b) high opacity.

Model complexity in terms of number of trainable parameters and the required tera floating point operations (TFLOPs) is shown in Table 6.

Table 6

Model complexity (UNet2D, UNet3D, and our).

6.3.

Lesion Quantification in the Internal Testing (N_cov = 197) and External Testing (N_ext = 98) Cohorts

In Table 7, we present the interquartile range (IQR) and coefficient of determination ($R^2$) on volumes between expert and automatic segmentation along with the overall per-patient mean Dice score for both internal as well as external test datasets. In the internal test set ($N_{\mathrm{cov}}=197$), no significant difference between volumes of expert and automatic segmentations was observed for GGOs ($p=0.3612$). Similarly, no significant difference between volumes of expert and automatic segmentations was observed for GGOs ($p=0.1563$) or high opacities in the external test set ($N_{\mathrm{ext}}=98$).

Table 7

Our model performance on Ncov=197 and Next=68.

The Bland–Altman analysis on the internal test set demonstrated a low bias of 0.56 (Fig. 6) and 18.61 ml (Fig. 7) for GGOs and high opacities, respectively. Similarly, the Bland-Altman analysis on the external test set demonstrated a low bias of 7.16 (Fig. 8) and 2.92 ml (Fig. 9) for GGOs and high opacities, respectively. After further analyzing the anomalies (cases outside the limit of agreement) in the Bland–Altman plots (Figs. 6 and 7), we observed that some input scans were corrupted due to various reasons including motion artefacts, errors in expert annotations, etc., as shown in Fig. 10. Thus, a significant ($p<0.001$) difference in high opacity volumes between expert and automatic segmentations was observed in the internal test set.

Fig. 6

Expert and automatic quantification of GGO in testing cohort ($N_{\mathrm{cov}}=197$). (a) Bland–Altman plot and (b) best fitting regression line.

Fig. 7

Expert and automatic quantification of high opacity in testing cohort ($N_{\mathrm{cov}}=197$). (a) Bland–Altman plot and (b) best fitting regression line.

Fig. 8

Expert and automatic quantification of GGO in external unseen testing cohort ($N_{\mathrm{ext}}=68$). (a) Bland–Altman plot and (b) best fitting regression line.

Fig. 9

Expert and automatic quantification of high opacity in external unseen testing cohort ($N_{\mathrm{ext}}=68$). (a) Bland–Altman plot and (b) best fitting regression line.

Fig. 10

Samples of extreme outliers from Bland–Altman plots. Highlighted red rectangle and ellipse are the areas of mis-classifications. Blue indicates GGOs, and yellow indicates high opacities.

The internal testing cohort consisted of 30 contrast enhanced and 167 non-contrast CT scans. We observed no significant difference ($p=0.2137$) between the mean Dice scores calculated for segmentations from non-contrast and contrast-enhanced CT scans, which were ($0.8939\pm 0.0663$) and ($0.8801\pm 0.0682$), respectively.

6.4.

Diagnostic Comparison

We trained the same convLSTM model with and without additional $N_{\text{control}}^{\text{train}}=197$ controls and tested them with five-fold cross-validation. We also tested an additional $N_{\text{control}}^{\text{test}}=498$ unseen controls with the best performing model from five-fold cross-validation. The AUROC with and without NLST in training was 0.965 and 0.959, respectively, but they did not reach significance. However, McNemar’s test results (Table 8) show that the model trained with an additional $N_{\text{control}}^{\text{train}}=197$ NLST cases significantly increased the specificity at 95% sensitivity of the model. Thus, adding NLST controls to the training decreased the false positive rate in diagnosis. The overall per-patient mean Dice score also improved drastically, as shown in Table 8.

Table 8

Diagnostic performance on Ntotal=Ncov+Ncontrol=892 NLST patients.