Deep learning for predicting COVID-19 malignant progression

Dataset statistics

All data enrolled in this retrospective study is obtained from two hospitals in Wuhan, including Wuhan Pulmonary Hospital and three branches of Tongji Hospital. The collection and use of these data from two hospitals were approved by their IRBs, respectively. 1,040 patients with mild COVID-19 pneumonia at admission are considered in our study, including 491 males and 549 females, aged 18 to 95 (57.51 ± 14.75). 32.3% of patients (336/1,040) malignantly progressed to a severe/critical stage during the hospitalization, while the remaining 67.7% (704/1,040) did not (Fig. 1). The selected data is divided into three cohorts, of which the cohort one is used for the single-center study, and the cohort two and the cohort three are used for the multicenter study. The clinical data in three cohorts is summarized in Table 1.

• Download figure
• Open in new tab

Fig. 1 Flowchart of patient selection

A total of 1,040 out of 2,742 patients are selected according to the inclusion criteria. All the 1,040 patients have the complete clinical data required for the study and 57.9% of them underwent serial chest CT imaging. RR: respiratory rate, SpO₂: blood oxygen saturation, PaO₂: arterial oxygen partial pressure, FiO₂: fraction of inspiration oxygen.

Performance evaluation and results

Our work advocates the use of a sequence of CT scans, captured at different timings after hospitalization, for accurate malignant progression. Unlike^8,9, we do not quantify CT scans to avoid information loss, but use a deep learning model to process the raw data directly. In the meantime, an effective integration of CT scans and clinical information underpins our system. The pipeline of our system is given in Fig. 2.

• Download figure
• Open in new tab

Fig. 2 The pipeline of our system about the prediction of COVID-19 malignant progression

First, 3D ResNet and MLP encode chest CT scans and clinical data, respectively. Then, we combine the two features and feed them into an LSTM to model the temporal information. Finally, several fully connected layers are exploited to make the prediction. CT: computed tomography, LSTM: long short-term memory, MLP: multilayer perceptron.

In Table 2, we compare the performance of our model against different methods, including linear discriminant analysis (LDA)¹³, support vector machine (SVM)¹⁴, and multilayer perceptron (MLP)¹². In this experiment, we adopt patients with mild symptoms and without malignant progression as the negative reference standard and divide the cohort one into the training cohort (80%) and the validation cohort (20%) randomly. We use five-fold cross-validation to evaluate our model. Our system, which fuses clinical information with the sequential CT scans, achieves a mean AUC of 0.920 (95% Cl: [0.861, 0.979]) and outperforms the best traditional machine learning methods (mean AUC of 0.767, 95% Cl: [0.725, 0.799]) by a large margin.

According to the results in Table 2, we can draw the following conclusions: 1) CT scans turns out to be more effective than quantitative CT features. For instance, with LSTM, using CT scans obtain a relative improvement of 2.3% over using quantitative CT features in AUC (0.920 vs. 0.899). Without LSTM, the improvement is more distinctive, reaching 8.1% (0.851 vs. 0.787). 2) Modeling the temporal information brings measurable benefits for boosting the performance of our system. This has been evidenced by an improvement of 14.2% (0.899 vs. 0.787) via using LSTM when quantitative CT features are used. Meanwhile, when CT scans are used, the improvement is 8.1% (0.920 vs. 0.851) with a difference of 0.069 in AUC. The corresponding ROCs in Fig. 3a further support our method.

• Download figure
• Open in new tab

Fig. 3 Comparison of ROC curves among different methods on three cohorts

Numbers before brackets are AUCs. Numbers in brackets are confidence intervals. a, ROC curves comparison on cohort one. LDA: linear discriminant analysis, SVM: support vector machine, MLP: multilayer perceptron, LSTM: long short-term memory. b, ROC curve on cohort two. PT: pre-trained. c, ROC curves comparison on cohort three. PT: pre-trained, FT: fine-tuning, DA: domain adaptation. In each class, the same number of samples are used during the domain adaptation process.

Performance evaluation in the multicenter study

A high-quality labeling process typically requires time-consuming human effort, which is a prominent drawback during the outbreak of COVID-19 where fast analysis is essential. Hence, how to use a small amount of labeled data to improve the generalization power of the system is of great practical significance. Inspired by metric-learning based methods¹⁵, we propose a domain adaptation method to adapt our model to a new domain with only a few labeled samples available. The details are given in the Methods section.

As Table 3 shows, when directly evaluating the model trained from the cohort one on the cohort two, a satisfactory performance is achieved (AUC: 0.885, 95% CI: [0.847, 0.923]). This is because the two cohorts are from different branches of the same hospital, which means that their data distributions are similar to some degree. However, when directly evaluating the model trained from the cohort one on the cohort three, the performance drops a lot. A mean AUC of 0.651 (95% CI: [0.558, 0.745]) is achieved, because the two cohorts are from different hospitals. When 10 labeled samples in the target domain are used, directly finetuning the model achieves a mean AUC of 0.738 (95% CI: [0.655, 0.820]), still inferior to our system (AUC: 0.862, 95% CI: [0.789, 0.935]). The corresponding ROCs in Fig. 3b,c further support our domain adaptation method.

Prognostic factors of clinical data

We further investigate the clinical indicators that contribute to the prediction of the malignant progression by adding a self-attention layer before the first fully connected layer of the MLP. As shown in Supplementary Figure 2, it automatically learns the attention weight corresponding to each clinical indicator, and each attention weight is normalized to (0, 1) by a sigmoid function to measure the importance of the related clinical indicator for the prediction task. The top 20 clinical indicators with the highest attention weights are listed in Supplementary Figure 3. The most important prognostic clinical indicators are myocardial injury (Troponin and Brain natriuretic peptide), followed by hepatic injury (Aspartate aminotransferase, Albumin, and y-Glutamyl transpeptidase), renal failure (Creatinine), and inflammatory status (Hypersensitive C-reactive protein, White cell count, CD3+CD4+T cells count, fever).

Clinical data acquisition and preprocessing

In Wuhan Pulmonary Hospital, data from 199 patients were collected from January 3, 2020 to February 13, 2020. In Tongji Hospital, data from 2,543 patients were collected from January 13, 2020 to March 16, 2020, of which 544 patients came from the Zhongfa branch, 363 patients came from the Guanggu branch, and 71 patients came from the Main branch. All patients were confirmed by a positive viral nucleic acid test. A subset of 1,040 adult patients belonging to the mild type at admission assessments is selected for further investigation. The inclusion criteria are all of the followings: 1) respiratory rate < 30 breaths per min; 2) resting blood oxygen saturation > 93%; 3) the ratio of arterial oxygen partial pressure to fraction of inspiration oxygen > 300 mm Hg; 4) non-ICU patients without shock, respiratory failure, mechanical ventilation, and failure of other organs. Anyone who fails to fulfill one of the criteria is considered to progress into a severe/critical stage according to the guidelines for the diagnosis and treatment of COVID-19 infection by the Chinese National Health Commission (Version 7)³⁰. The medical history, results of physical examination, and laboratory tests were all collected from the HIS system. The time points of symptoms onset and the beginning of severe/critical stage are recorded for further selections of available CT scans. All the patients in the Main branch have mild prognoses, so our research excludes patients from this branch. Furthermore, considering the age imbalance problem, we exclude patients under 18 years old. We finally obtain 61 clinical indicators for each sample.

CT data acquisition and preprocessing

All the patients underwent serial pulmonary CT exams on dedicated CT scanners (GE, SIEMENS, TOSHIBA, and UNITED IMAGING) in two hospitals with the following parameters: slice thickness 1-3 mm, slice gap 0 mm, 130 kV, 50 mAs. All the CT scans before the severe/critical stage are included to segment the masks of bilateral lungs and pneumonia on an autonomous system (HUAWEI CLOUD Launches AI-Assisted Diagnosis Platform for COVID-19). Each CT scan is downsampled to a 64 × 64 × 64 (width x height x slice) tensor. CT image values are clamped to the range [-1250, 250]. Data augmentations including random horizontal flips and random rotations are used.

As traditional machine learning methods cannot handle raw CT scans directly, we design a set of hand-crafted quantitative features for experimental comparison. These quantitative features include the infection pixels proportion and the average CT value of infection pixels for each CT slice. Finally, a 128-dimensional vector can be obtained for each CT scan. We use zero-padding for the missing CT scan tensor.

Network architecture and training process

Fig. 2 illustrates the pipeline of our system. As is shown, the input of our model includes the clinical data and a sequence of CT scans obtained at different time points. Specifically, the clinical data is a 61-dimensional vector which is processed by an MLP with identity connections (Supplementary Figure 5). In addition, each CT scan is encoded into a 128-dimensional feature vector by 3D ResNet.

To model the temporal information across the sequence of CT features, we use LSTM for its high capacity in modeling such information³¹ and densely combine the clinical feature and the CT feature via concatenation at each time step. LSTM employed in this study is a single-layer network with an embedding dimension of 189 and a hidden dimension of 378. The output of LSTM, a 378 × 7 tensor, is flattened and then fed into several fully connected layers. Finally, we normalize the output with a softmax layer, which can be interpreted as the probability of the patient’s conversion to the severe/critical stage. The whole model is trained with the cross-entropy loss. The detailed architecture of our model is given in Supplementary Table 1.

Domain adaptation process

In domain adaptation, we use a metric-based method by using a few labeled samples to bridge the domain gap between the source center (cohort one) and the target center (cohort three). Specifically, our method can be decomposed into two stages: a pre-training stage and a domain adaptation stage. For the pre-training stage, we first train a model on the source center then remove the classifier to get the pre-trained encoder f_ϕ. For the domain adaptation stage, which is the core of our method, we adapt the pre-trained model through a metric-based approach, passing the prototype representations learned from the source center to the target center. The details of this stage are elaborated as follows.

First, we randomly select N labeled samples from each class in the target center as the support-set to compute prototypes. At the same time, we randomly choose one sample per class in the target center as the query-set to compute distances to the prototypes in the embedding space. Specifically, in each class let S = {(x₁, y₁), …, (x_N,y_N)} denotes a small support-set of N labeled samples, where x_i is the feature vector of an example and y_i ∊ {0, 1} is the corresponding label. S_k denotes the support-set of examples with class k ∊ {0, 1}. We compute the mean vector w_k of the embedded support points as prototypes for the two classes:

Second, we compute the predicted probability distribution for each sample in the query-set based on a softmax over its cosine similarities (denoted by sim(.)) with the prototypes in the embedding space:

Similar to these studies^32-34, we also use a learnable parameter τ to control the probability distribution sharpness generated by the softmax function during training.

Finally, we train our framework with the classification and similarity losses jointly. Concretely, the classification loss is computed based on p and the labels in the query-set: where L_CE is the cross-entropy loss. The similarity loss is adopted to increase the distance between the two class prototypes, which is defined as:

Based on the above two losses, our final loss function is formulated as: where λ is the coefficient to balance the two loss terms.

In the testing stage, all labeled samples in the support-set and the query-set are used to compute prototypes, a test sample will be classified by the similarity between the feature vector of the sample and the prototype of class k.

Implementation details

The proposed network is implemented using Python (Version 3.6 with scipy, scikit-learn, and PyTorch). For single-center experiments on the cohort one, the network is trained by Adam optimizer with an initial learning rate of 0.05, and a batch size of 32 on a single NVIDIA Titan X GPU. The learning rate is decayed by a factor of 10 every 30 epochs. We train our model for 100 epochs. Model weights are initialized with Kaiming method³⁵ and bias are initialized as 0. For multicenter experiments on the cohort three, we use the same optimizer settings as single-center experiments, but with a fixed learning rate of 0.01. We finetune 50 epochs for domain adaptation and the batch size is the same as the number of labeled samples with a maximum of 20. The cosine scaling parameter τ is initialized to 5 with a fixed learning rate of 0.1. The loss coefficient λ is 0.5. The margin of the similarity loss is set to 0.2.

Performance evaluation criterion

The AUC, accuracy, sensitivity, specificity, and ROC curve are used to evaluate the model performance. The calculation method is shown in Supplementary methods. The 95% bilateral confidence interval is used for all metrics, where the AUC metric uses the Wald-cc interval^36,37 and the other metrics use the Wilson interval³⁸.