An artificial neural network classification method employing longitudinally monitored immune biomarkers to predict the clinical outcome of critically ill COVID-19 patients

Train/test dataset

One clinical cohort of COVID-19 ICU patients was considered for training/testing the model. This study took place in Cleveland, USA. Patients were enrolled from March 2020 to May 2020. Patients were hospitalized either because they were potential candidates for mechanical ventilation and/or because they were judged to be in an unstable condition requiring intensive medical or nursing care. This cohort was composed of 45 patients. The patients had longitudinal samplings on their biomarkers consisting of multiple time points (first day of admission, fifth day, eighth day, eleventh day, and 15th day) unevenly spread. All the samples associated with each patient were averaged and associated with a clinical outcome, i.e., deceased or recovered. A panel of biomarkers that flag a potential dysregulated immune response was systematically put together (Material S1) and profiled using the Ella Simple Plex Immunoassay (San Jose, CA, USA). The biomarkers used were as follows: Intercellular Adhesion Molecule 1 (ICAM-1), Lipocalin-2 (LIPO), Myeloperoxidase (MPO), Vascular Cell Adhesion Molecule 1 (VCAM-1), D-Dimer, E-selectin (E-SEL), Ferritin, Surfactant Protein (SP-D), Programmed Death-Ligand 1 (PDL1), Granulocyte Colony-Stimulating Factor (G-CSF), Interleukin 1 beta (IL-1b), Vascular Endothelial Growth Factor C (VEGFC), Angiopoietin 2 (ANG2), C − X − C Motif Chemokine Ligand 10 (CXCL10), Granulocyte Macrophage Colony-Stimulating Factor (GM-CSF), Interleukin 10 (IL-10), Interleukin 17 (IL-17A), Interleukin 1 receptor antagonist (IL-1ra), Interleukin 6 (IL-6), Interleukin 7 (IL-7), C-C Chemokine Ligand 2 (CCL2), Granzyme B (GRANB), Interferon gamma (IFNg), Interleukin 12 (IL-12), Interleukin 15 (IL-15), Interleukin 2 (IL-2), Interleukin 4 (IL-4), and Tumor Necrosis Factor-alpha (TNF- α). The train/test dataset is available in Material S2.

Validation dataset

A second clinical cohort of ICU COVID-19 patients was considered for validating the model. This cohort took place in Dublin, Ireland. Patients (n = 50) were recruited from September 2020 to March 2021. The patients were binarized to 30 recovered and 20 deceased. We have opted to use these patients to validate the model because all patients were critically admitted to the ICU. This cohort was composed of the same 28 biomarkers of the training/testing cohort that were quantified through the same Simple Plex Immunoassay (San Jose, CA, USA). The patients were sampled during their ICU admission (1st day), day three, and 14 days later. The concentration of the available samples was averaged so each patient had one unique value per biomarker. The validation dataset is available in Material S3.

Dimensionality reduction

The two clinical cohorts employed in this study have concentration levels of 28 biomarkers. It has been reported that a dataset with many dimensions will eventually decrease the performance of a given machine learning algorithm or might produce biased results as the fitting curve is too complex (Taylor, 2019). To reduce the number of variables to characterize a patient, a principal component analysis (PCA) was achieved through the prcomp function implemented in the R stats package (version 4.1.2). The loading scores of the principal components 1 and 2 were obtained in the $ rotation element that a prcomp object has. The R code (R version 4.1.2) used for reducing the dimensions of the cohorts of this study is found at http://dx.doi.org/10.5281/zenodo.6643238.

Classification

We have benchmarked five algorithms for classifying the patients from the train/test and validation datasets according to their labels (i.e., 0 deceased and 1 recovered, included in Material S2 and Material S3): support vector machines (SVMs), random forest (RF), classification and regression trees (CART), k-nearest neighbors (KNN), XGBoost, and ANNs. The feasibility of the classification of each algorithm was measured through the accuracy metric upon a 5-fold cross validation dataset. The SVM, RF, CART, and KNN were built in R and the script is available at https://github.com/gustavsganzerla/covid-biomarker/blob/main/different_classif-models.R. The XGBoost classifier was built in Python with the XGBClassifier library (version 1.4.0) and is available at https://github.com/gustavsganzerla/covid-biomarker/blob/main/xgboost_biomarkers_COVID.ipynb.

A deep learning instance of ANNs was used. The ANN simulations took place in the Python language through the TensorFlow library in its version 2.8.0. To scale the entry data under the same magnitude, a normalization process was conducted through the StandardScaler function found in the sklearn.preprocessing library. To validate the simulation process, a cross-validation process was considered, where k =5. This method involves reserving 4/5 of the dataset for the training process and the remaining 1/5 for testing. Then, the process is repeated k times so every data point is covered both in training and testing. The training data consists of 80% of the data.

The ANNs were trained, tested, and validated (the latter with external data to grant generalization capacity to the model). To assess a binary prediction, a confusion matrix was built. Deceased patients have the 0 label while recovered patients have 1. This representation of the classification assigns a prediction to a: (i) True Positive (TP), i.e., patients who survived and the model classified as survivors, (ii) True Negative (TN), i.e., patients who deceased and the model indicated so, (iii) False Positive (FP), patients who actually died and the model indicated they survived, and iv) False Negative (FN), i.e., patients who survived and the model indicated they died. The assessment of the model involves calculating the following performance metrics.

Accuracy, which measures the number of correct predictions in the whole dataset; its calculation is achieved by Eq. (1). (1)${\displaystyle Accuracy=\frac{\left(TP+TN\right)}{\left(TP+TN+FP+FN\right)}.}$ Precision, which measures the ratio between the TPs among all the positives, its calculation is achieved by Eq. (2). (2)${\displaystyle Precision=\frac{TP}{\left(TP+FP\right)}.}$

Recall, which measures the percentage of the model identifying TPs; Eq. (3) calculates it. (3)${\displaystyle Recall=\frac{TP}{\left(TP+FN\right)}.}$

Finally, the specificity metric calculates the detection rate of TNs throughout the entire dataset. It is obtained through Eq. (4). (4)${\displaystyle Specificity=\frac{TN}{\left(TN+FP\right)}}$

Next, the model had its recall and specificity scores assessed to obtain the Receiver Operator Characteristic (ROC) curves, which measured the trade-off between recall and specificity following different thresholds applied to the outcome of the sigmoidal function used in the output neuron. The default threshold value (0.5) was preserved in all instances of the ANN simulation. A second instance of the classification procedure was achieved by isolating the sex of the patients. The Python code that implemented all the ANNs simulations in this study is available at https://github.com/gustavsganzerla/covid-biomarker.

Ethics approval

The training/testing cohort obtained an Institutional Review Board (IRB) Approval from MetroHealth Medical Center in Cleveland, Ohio IRB 20-00198 on March 25, 2020. The validation cohort obtained an IRB Approval from SJH/TUH Joint Research Ethics Committee and The Health Research Consent Declaration Committee (HRCDC) under the register REC: 2020-05 List 17 on March 2, 2020.

Dimensionality reduction

To reduce the dimensionality of the dataset and consequently have fewer neurons being used in the input layer, we applied a Principal Component Analysis (PCA). From it, we selected the first two principal components (PCs) (Fig. 1A); then we measured the loading score in these two components and extracted the variables that most contributed to explain the components’ variance (Fig. 1B). PC1 is responsible for 22.7% of the data variance while PC2 is responsible for15%. We maintained this proportion and selected three biomarkers from PC1 and two from PC2. Therefore, we report a reduced number of biomarkers that explains the variability found in the dataset; i.e., VCAM.1, TNF- α, ICAM.1, CXCL10, G.CSF. In addition, we have also screened the 5 worst-ranking biomarkers by the PCA for validation purposes (we maintained the same proportion as that used in the 5 best-ranking PCA biomarkers).

Defining the classification procedure

We selected RF, SVMs, CART, XGBoost, KNN, and ANNs to classify the data. As the accuracy performance of the first four techniques did not show satisfactory results in correctly labeling the train/test and validation cohorts (Material S4), we therefore opted to use Deep Learning ANNs as the previous models did not show a capacity of generalizing upon unseen data.

To determine the optimal architecture of ANNs to be used in this work, a series of parameters were set. First, the input layer of the ANN contains the number of biomarkers that will be employed to classify patients’ outcomes and nature of infection. Next, we created different ANN models to determine the optimal number of hidden layers and their number of neurons (Material S5). From that, we opted to use an architecture consisting of four hidden layers with 20 neurons each, where the neurons are activated by a ReLU function. The loss function was set to binary cross entropy and the Adam optimizer was chosen. We allowed the weights of the ANN to be updated 200 epochs; this was found to be the limit in which the error did not drop any further. Since all the further prediction steps deal with binary classification, we set one neuron in the output layer; this neuron is activated by a sigmoid function since it yields a probability (i.e., probability of recovering or dying).

Classifying patients’ outcome with 28 biomarkers

We predicted the clinical outcome of COVID-19 patients using 28 biomarkers as the input of the ANN model. The results (Fig. 2A) suggest a satisfactory classification of the patients. We further explored our classification metrics by providing the ROC curve (Fig. 2B) in which we adjusted the decision threshold obtained as the output by the sigmoidal function in the output neuron. Lastly, we show in Fig. 2C, the error dropping with increasing epochs.

An optimal set of biomarkers succeeds in classifying patients’ outcomes

To test the prediction capacity of the biomarkers identified by the PCA, we ran an ANN with five neurons in the input layer. To validate the prediction power of the ANN with the PCA biomarkers, we performed two classification procedures with different input sets: one with the five best-ranking PCA biomarkers (three from PC1 and two from PC2) in the input layer and a second with the five worst-ranking PCA biomarkers (three from PC1 and two from PC2). In Fig. 3, we show the distinctive classification obtained by both sets of biomarkers (Fig. 3A). Next, we show the ROC curves for each point in the decision threshold (Fig. 3B). Finally, in Fig. 3C, we show the error drop rate, during the execution of both ANNs. We found the ANN trained with the 5-best ranking PCA biomarkers ANN dropped its error at a significantly different rate than its counterpart (p = 2.2e−16) with the increase of epochs. In addition, due to the time of enrollment, the patients from the train/test cohort were likely infected with the Wuhan D614G strain of the virus.

We combined male and female patients of the two cohorts in separate datasets to assess the classification capacity gender-wise (Table 1). We maintained the same input biomarkers identified by the PCA analysis. We found that the classification of female patients outperforms the metrics of male patients despite the accuracy of male patients being higher than that of females (i.e., 77% vs. 73.66%, respectively).

The prediction capacity of individual biomarkers

To test the clinical outcome prediction capacity of each biomarker, we ran individual ANNs with a single neuron as input. In these regards, we observed the ROC curves for each individual biomarker (Fig. 4) and compared it to how the error dropped after 200 epochs in the ANN model (Fig. 5). The same biomarkers with satisfactory ROC scores matched the low error after 200 epochs (Table 2).

The biomarkers that individually showed a satisfactory prediction capacity (balanced performance metrics) are TNF- α, PDL1, LIPO, ICAM1, and VCAM; from these, TNF- α, ICAM.1, and VCAM are the main contributors of the PC1 found by the PCA analysis.

Validation of the findings with an external cohort

In order to provide validation to the models previously identified, an external source of COVID-19 patients. At that time, the circling variant in the British Isles was B.1.1.7. In Fig. 6A, we assess the generalization capacity of three different classification approaches explored in this work. From it, the patients were divided into two classes (i.e., recovered and deceased), and each class was validated with the classification models. Firstly, the classification with 28 biomarkers only classified the deceased patients correctly. In the second approach, the five best ranking PCA biomarkers produced a balanced prediction of the two classes. Finally, the classification with the best individual biomarker (i.e., TNF- α) also only managed to yield a satisfactory classification of deceased patients. In Fig. 6B, we isolate the model that well generalized the external data (i.e., five best-ranking PCA biomarkers) and map the validation data into the ANN’s sigmoid function for which both deceased and recovered patients were satisfactorily classified. The average function outcome is 0.42 and 0.61 for deceased and recovered, respectively.