Identifying Immunological and Clinical Predictors of COVID-19 Severity and Sequelae by Mathematical Modeling

Since its emergence as a pandemic in March 2020, coronavirus disease (COVID-19) outcome has been explored via several predictive models, using specific clinical or biochemical parameters. In the current study, we developed an integrative non-linear predictive model of COVID-19 outcome, using clinical, biochemical, immunological, and radiological data of patients with different disease severities. Initially, the immunological signature of the disease was investigated through transcriptomics analysis of nasopharyngeal swab samples of patients with different COVID-19 severity versus control subjects (exploratory cohort, n=61), identifying significant differential expression of several cytokines. Accordingly, 24 cytokines were validated using a multiplex assay in the serum of COVID-19 patients and control subjects (validation cohort, n=77). Predictors of severity were Interleukin (IL)-10, Programmed Death-Ligand-1 (PDL-1), Tumor necrosis factors-α, absolute neutrophil count, C-reactive protein, lactate dehydrogenase, blood urea nitrogen, and ferritin; with high predictive efficacy (AUC=0.93 and 0.98 using ROC analysis of the predictive capacity of cytokines and biochemical markers, respectively). Increased IL-6 and granzyme B were found to predict liver injury in COVID-19 patients, whereas interferon-gamma (IFN-γ), IL-1 receptor-a (IL-1Ra) and PD-L1 were predictors of remarkable radiological findings. The model revealed consistent elevation of IL-15 and IL-10 in severe cases. Combining basic biochemical and radiological investigations with a limited number of curated cytokines will likely attain accurate predictive value in COVID-19. The model-derived cytokines highlight critical pathways in the pathophysiology of the COVID-19 with insight towards potential therapeutic targets. Our modeling methodology can be implemented using new datasets to identify key players and predict outcomes in new variants of COVID-19.

1 Introduction

The coronavirus disease (COVID-19) has been following a non-linear evolution through the pandemic, starting with one variant that mutated into at least four dominant subtypes. Early prediction of COVID-19 outcome is crucial to direct resource allocation by the health care system and to triage the patients to receive the optimum clinical management. Despite the broad spectrum of presentations, a significant turning point in the course of the disease is the development of abrupt systemic elevation of a myriad of inflammatory cytokines and chemokines (the cytokine storm- CS). In this phase of COVID-19, multiple organ failure progressing to circulatory shock is the leading cause of death. The CS is accompanied by a myriad of biochemical and radiological findings (1). A key determinant factor of COVID19 progression is the uncontrolled dysregulation production of cytokines and chemokines, resulting in the development of a cytokine storm, systemic inflammation, and consequently multi-organ failure (2). The presence of the cytokine storm was associated with COVID-19 severity as previously reported (3, 4), where the serum levels of cytokines in COVID-19 patients were significantly correlated with the severity of the disease and acted as warning indicators of the severity and progression of COVID-19.

Interpreting the role of cytokines, their predictive value and therapeutic potential is still a significant challenge in the context of COVID-19. An example of an incomplete understanding of CS and its pathogenesis is the uprise and drop of tocilizumab. As interleukin-6 (IL-6) is a critical cytokine in CS-induced mortality in patients receiving engineered T cell therapy, it was first suggested as a potential therapeutic target for COVID-19 CS. However, a randomized, double-blind Phase III COVACTA trial failed to reveal a significant reduction in mortality by using tocilizumab in COVID-19 (NCT04320615) (5), mandating further (re)search for additional key players in the CS pathogenesis.

The COVID-19 is an adaptive dynamic disease that has witnessed the SARS-CoV-2 mutated multiple times since March 2020. It is highly expected that SARS-CoV-2 will persist as an endemic infection, with epidemic peaks (6), as witnessed with the 4th and 5th waves in some countries. Many of the early models for COVID-19 failed to predict many aspects of the disease (7). Part of the issue is that COVID-19 is a non-linear disease. Many molecular studies were carried out to understand COVID-19 initiation and progression. However, such studies faced various challenges, including the curse of dimensionality (where the total number of severely infected patients is relatively low but each patient has a high number of data points) and inability to find optimal solutions across the general problem and thus end up with sub-solutions (local minima) (8). Artificial Intelligence (AI) is designed to find global solutions to multi-dimensional data. In the context of COVID-19, AI offers vital tools to find better predictors. However, AI has a few limitations in biomedical applications, mainly because AI solutions can be skewed by noise and thus requires well-annotated datasets with a clear understanding of the measured parameters. Integrating clinical, radiological and biochemical tests is highly recommended to achieve the ultimate benefit of modeling the disease.

Interestingly, stochastic modeling was previously used to model the human immune response to the yellow fever vaccine (9). Since COVID-19 is linked to immune response, modeling of the SARS-CoV-2 infection have been extensively published on different aspects of the disease, including the immune system using multiple ODEs to model immune cells, antibodies and cytokines (10–13), and on the clinical and radiological data (14–16). A few models on cytokine release syndrome in other diseases were also created (17–19). Investigating the immune response signature in COVID-19 yielded various biomarkers in different studies. A previous retrospective analysis suggested IL-6, IL-8 and TNF-α as independent predictors of patient survival (20). More recently, Perreau et al., 2021 suggested hepatic growth factor (HGF) and CXCL13 as predictors of severity and mortality of COVID-19 (21).

In the current study, we hypothesized that integrative analysis of pro-inflammatory, anti-inflammatory cytokines and checkpoint markers in addition to key clinical, biochemical, and radiological parameters could predict COVID-19 outcomes with higher predictive accuracy than individual parameters. Guided by the transcriptomics analysis of nasopharyngeal swabs, we curated a panel of 24 cytokines to be assayed using a multiplex assay with high intra- and inter-assay precision to reflect the immune response in our model, using a small amount of serum. Added to the 24 entries of cytokine levels, we also included 63 entries of clinical, biochemical and radiological parameters of well-characterized patients. In this study, we are introducing our clinically applicable integrative model as a predictive tool for COVID-19 severity and sequelae that will hopefully help guide clinical decision and management strategies. In addition, the study highlights potential therapeutic targets via identifying key players in the cytokine storm.

2 Patients and Methods

2.1 COVID-19 Patients’ and Healthy Controls’ Criteria

Nasal swab samples were collected from 50 COVID-19 patients (10 Asymptomatic, 11 mild, 13 moderate, and 16 severe patients; SARS-CoV-2 infections is confirmed by PCR), in addition to 11 healthy donors, at Rashid Hospital in Dubai, following the approval of the ethical committee at Dubai Health Authority (DSREC-04/2020_09). All patients were recruited between February-March 2020, and hence did not receive COVID-19 vaccine.

Peripheral venous blood samples of 37 COVID-19 patients were collected following the approval of the ethical committee at Rashid Hospital in Dubai (DSREC-04/2020_19), in addition to 40 healthy controls. All the included patients were recruited between June-July 2020, and hence did not receive COVID-19 vaccine.

Patients were classified into the respective group severity as follows: (A) Mild-moderate: no or mild pneumonia, (B) severe: patients with at least one of the following symptoms: shortness of breath (breathing rate ≥ 30/min), SaO2 at rest ≤ 93%, partial pressure of oxygen in arterial blood (PaO2)/inspired oxygen fraction (FiO2) ≤ 300 mmHg, or lung infiltrates > 50% within 24 to 48 h. Clinical and biochemical data were collected. Also, computed tomography (CT) imaging was performed, followed by an assessment using the COVID-19 Reporting and Data System (CO-RADS) as a standardized assessment of pulmonary involvement of COVID-19 (22). The Co-RAD categories correspond to the corresponding level of suspicion of pulmonary involvement in COVID-19. 0= scan is technically insufficient to assigning a score; 1=normal or non-infectious; 2= typical for other infection but not COVID-19; 3 = features are compatible with COVID-19 but also other diseases; 4 = highly suspicious for COVID-19; 5 = typical for COVID-19; and 6 = RT=PCR positive for SARS-COV-2.

The healthy controls (age: 47.18 ± 16.6 years, 24 males and 16 females, BMI: 25.9 ± 3.11 Kg/m²) were filtered from an initial cohort of 150 controls to include only those with a non-obese BMI and normal HbA1c to avoid having any confounding factors such as obesity or prediabetes.

2.2 Whole Transcriptome and Bioinformatics Analysis of Nasal Swab Samples From COVID-19 Patients and Healthy Controls

One ng of RNA of each sample was analyzed using targeted whole RNA-seq with AmpliSeq whole transcriptome on S5 system (Thermo Fisher Scientific), and RNA-seq data were analyzed using the Ion Torrent Software Suite version 5.4. Alignment was carried out using the Torrent Mapping Alignment Program (TMAP), as described in (23).

The expression counts of the nasal swap samples of COVID-19 patients and healthy controls were normalized across samples using the DESeq2 normalization method. Differentially expressed genes between each of the severity groups against the healthy control group were identified using the Bioconductor package DESeq2. Differentially expressed genes with adjusted p-value <0.05 and fold change >2 or <0.5 were considered statistically significant. The adjusted p-value was calculated using false discovery rate (FDR) according to Benjamini Hochberg method (24).

2.3 Bioinformatics Analysis of Publicly Available COVID-19 Whole Blood RNA-Seq Dataset

In addition, whole blood bulk RNA sequencing data (Normalized counts) deposited by Bernardes et al. (25) were retrieved from: https://github.com/Systems-Immunology-IKMB/COVIDOMICs/tree/main/TF_enrichment/TF_enrichment_analysis-main/data. The dataset included samples from 42 COVID-19 patients (12 asymptomatic, 11 mild, 6 complicated, 4 complicated incremental, 6 complicated hyper-inflammatory, and 3 critical patients), in addition to 14 healthy donors. Statistical significance of the differential expression of cytokines between the disease severity groups was analyzed using one-way ANOVA with post hoc Tukey’s multiple comparisons test in Graph Pad Prism (version 5.01). A p < 0.05 was considered statistically significant.

2.4 Collection of COVID-19 Patients’ Blood Samples

Peripheral Venous blood samples of 37 COVID-19 patients were collected following the approval of the ethical committee at Rashid Hospital in Dubai (DSREC-04/2020_19). Ethylene-diamine-tera-acetic Acid (EDTA) containing tubes were used to collect the blood samples. Then, serum was separated for the cytokine assays Forty blood samples were obtained from before the first case of COVID-19 infection in the UAE (MO-HAP/DXB/SUBC/No.14/2017).

2.5 Cytokine Assay

Given the results of previous steps, various cytokines of significance were assessed in the sera (50ul sample) of the COVID-19 patients and healthy controls using the Human Immunotherapy Magnetic Luminex Performance Assay 24-plex Fixed Panel (R&D systems, USA). The assessment was done using the Bioplex-200 system (Biorad, USA). A list of the curated cytokines is provided in Supplementary Table S1.

2.6 Statistical Analysis

Groups with different severity were compared after testing normality (Kolmogorov-Smirnov and Shapiro- Wilk tests). If the p-value is <0.05, non-parametric tests were used (Mann Whitney for comparing two groups, or Kruskal-Wallis to compare more than two groups). We grouped the mild and moderate as (non-severe). Of all variables, only age, BMI and platelet count followed normal distribution where the unpaired t-test was used. Power calculation was performed based on Wei et al. (26), setting the statistical power at 0.8, α = 0.05, and using the mean values of different cytokine levels in mild-moderate versus severe cases. The minimum number required in each group was estimated to be 15. The statistical package SPSS (v.28) was used for statistical analyses.

2.7 Machine Learning

Machine learning was used to reduce the set of clinical parameters and identify the optimal set of parameters to stratify the patients according to the different aspects of COVID-19 pathogenesis. A mixture of unsupervised hierarchical and K-means clustering analyses were performed in R (code in the Supplementary Material) to assess the separation of COVID-19 cases according to the blood protein expression levels of cytokine quantified using BIO PLEX-200. k = 6 was used as the number of clusters for the k-means clustering analysis. The k-means by storing all the labeled examples, and using them directly for inference on new data.

2.8 Mathematical Modeling

Mathematical modeling was used to identify key cytokines and biochemical markers that can stratify the clinical parameters collected in the study. The mathematical modeling was carried out using two models that were integrated subsequently. The first is ANOVA multivariate model with Bonferroni’s multiple testing. This is used to identify the variables that are significantly different amongst the various compared patients’ groups as well as different parameters denoting severity (e.g., mechanical ventilation, radiological findings, complications, e.g. liver injury).

The second model is the Stepwise linear regression model. This dynamic method systematically reduces the set of parameters (e.g., cytokines, biochemical parameters), depending on the significant interaction between the variables. The ANOVA multivariate model with Bonferroni’s multiple testing and the Stepwise linear regression model. Both models were applied to a combination of categorical (e.g., disease severity, oxygen support) and continuous data (e.g., protein expression and level of biochemical markers). Two-sided p <0.05 were considered to be statistically significant. ROC analysis was performed to assess the predictive efficacy of the predictors identified from the two mathematical models.

Different parameters were measured to check the model accuracy. if a is true positive, b is false positive, c is false negative, and d is true negative, the sensitivity was calculated as [a/(a+c)]×100; specificity as [d/(b+d)]×100; positive predictive value as [a/(a+b)]×100; and negative predictive value as [d/(c+d)]×100. Positive likelihood ratio was calculated as Sensitivity/(1-Specificity); negative likelihood ratio as (1- Sensitivity)/Specificity (26, 27).

3 Results

3.1 Nasopharyngeal Samples Identify Cytokines as Top Upregulated DEGs and Signaling Through Cytokines as Top Upregulated Pathway

Previous studies on nasopharyngeal swab samples were highly insightful on shifts in the immune landscape in association with COVID-19 (28), in contrast to the transcriptomic signature associated with different types of respiratory infections (29). However, general shifts in transcriptomic profiles associated with COVID-19 severity warranted further dissection and biological validation. Therefore, we carried out a transcriptomics analysis of nasopharyngeal samples collected from asymptomatic, mild, moderate, and severe patients using samples from healthy donors as a reference (Table 1). The upregulated transcriptome was significantly enriched in cytokine signaling and immune response pathways in moderate and severe COVID-19 patients (Figure 1 and Supplementary Figure S1), with several cytokines being in the top 100 DEGs. Our analysis revealed the significant upregulation of genes expressing IFN-γ, CXCL10, IL-33, Granzyme-B, and PD-L1 in moderate COVID-19 patients only and IL-8, IL-1Ra, IFN-α, CCL4, TNF, CCL3, and IL-1ß was in moderate and severe COVID-19 patients; in comparison, to healthy donors or asymptomatic patients (Figure 2A).

TABLE 1

Table 1 Demographic and Clinical data of the exploratory cohort (COVID-19 patients tested by transcriptomics analysis of nasopharyngeal swabs).

FIGURE 1

Figure 1 Pathways Enrichment is the nasopharyngeal swab samples of moderate and severe COVID-19 patients. Functional clustering and pathway analysis of the significantly upregulated genes in the nasopharyngeal swab samples collected from (A) moderate and (B) severe COVID-19 patients in comparison to healthy patients. DEGs were identified using DESeq2 algorithm; the genes were filtered according to adjusted p-value of <0.05 and fold change >2 or <0.5. The functional clustering analysis was performed using Metascape; p-value cut-off for pathways inclusion was <0.01.

FIGURE 2

Figure 2 Transcriptomics Analysis of nasopharyngeal swab samples and whole blood samples from COVID-19 patients. (A) Gene expression of cytokines and inflammatory mediators from the nasopharyngeal swap RNA-seq data compared across the different severity groups of COVID-19 cases (asymptomatic, mild, moderate, and severe) in reference to the non-COVID-19 control group. The data represented as log 2 normalized expression, where the normalized was performed using DESeq2 normalization approach across all the examined samples. (B) Gene expression of cytokines and inflammatory mediators from the whole blood RNA-seq dataset, compared across the different severity groups of COVID-19 cases (asymptomatic, mild, complicated, and critical) in reference to the non-COVID-19 control group. The data represented as log 2 normalized expression. * represents p-value < 0.05; ** represents p-value < 0.01; *** represents p-value < 0.001; analyzed using one-way ANOVA with post hoc Tukey’s multiple comparisons test.

To investigate whether the upregulation of these cytokines and inflammatory mediators is localized or systemic, the COVID-19 patients whole blood RNA-seq publicly available dataset (25) was analyzed. The analysis revealed the significant upregulation of IL-10 in mild, complicated and critical cases; IL-15 in mild and complicated cases; PD-L1 in complicated and critical cases; and IFN-γ in mild cases (Figure 2B).

3.2 Cross-Validation of Cytokines Using Bio-Plex

3.2.1 Recruited Patients

Based on previous findings and existing knowledge from previous publications, we examined the association of these cytokines with different clinical aspects of COVID-19 pathogenesis in a new cohort of 37 patients (on day 0-5 of admission, followed up for 4 weeks), and 40 age- and gender-matched healthy controls [Out of initial 150 control subjects, we selected 40, with a non-obese BMI (mean= 25.9 ± 3.11 Kg/m²) and normal HbA1c (<5.8%), to avoid having any confounding factors that may affect the cytokine levels, such as obesity or prediabetes]. Patients’ clinical data is provided in Table 2.

TABLE 2

Table 2 Demographic, Clinical and laboratory data of the validation cohort (COVID-19 patients tested for cytokine).

3.2.2 Cytokine Assay

In view of the transcriptomics analysis, the curated panel included cytokines (inflammatory and anti-inflammatory), chemokines and other immune-related molecules such as checkpoint markers, receptors and cytotoxic mediators (Supplementary Table S1). Out of the 24 investigated cytokines, 17 markers showed a differential pattern in COVID-19 patients compared to healthy controls (Figure 3).

FIGURE 3

Figure 3 Cytokine assessment in healthy control subjects (n =40), mild-moderate COVID-19 (n= 20) and severe COVID-19 (n=17) patients. (A) Inflammatory, (B) anti-inflammatory cytokines, (C) chemokines, and (D) checkpoint markers, receptors and cytotoxic mediators were assessed in mild-moderate and severe COVID-19 patients and their levels compared to healthy controls. Data is expressed as mean ± standard error of mean (SEM). *p<0.05, ** p<0.01, ***p<0.001, and **** p<0.0001.

As shown in Figure 3A, the levels of the pro-inflammatory cytokines GM-CSF, IL-6, IL-15, and IFN-α were higher in mild-moderate COVID-19, compared to healthy controls, with a further increase in severe COVID-19. On another note, levels of TNF-α and IL-17A were similarly more elevated in the mild-moderate and severe COVID-19 patients than in healthy controls. Also, IL-1β levels were found to be increased in mild-moderate COVID-19 patients that were restored in severe patients. While the component IL-12p70 showed a reduction in the serum levels of COVID-19 patients with a significant observed decrease in the severe patients’ group, previous studies reported no difference in the plasma levels of IL-12p70 (30). On the other hand, anti-inflammatory cytokines such as IL-1Ra and IL-10 showed a sequential increase in mild-moderate and severe COVID-19, while IL-4 showed a significant increase in mild-moderate COVID-19 (Figure 3B).

The chemokines MCP-1 (CCL2) IP-10 (CXCL10) incrementally increased levels in mild-moderate and severe cases (contributors to pulmonary pathogenesis). MIP1β (CCL4) increased equally in mild and severe cases. IL-8 (CXCL8) increased mild-moderate cases and decreased in severe cases (but still significantly higher than normal controls), (Figure 3C). As illustrated in Figure 3D, PD-L1 was found to be higher in severe COVID-19 than healthy controls or mild-moderate COVID-19. The transmembrane glycoprotein CD40 ligand and the cytotoxic molecule granzyme B showed a significant increase in mild-moderate COVID-19 patients compared to healthy controls. However, they were reduced in the severe patients (but still significantly higher than normal controls).

3.3 Exploration of Cytokines Expression Levels in Association With COVID-19 Disease Severity Using Machine Learning Techniques

The protein expression data of the cytokines and inflammatory mediators were further explored with machine learning approaches to identify the optimal set of parameters to stratify the patients according to different aspects of COVID-19 pathogenesis. Initially, unsupervised hierarchical and k-means clustering were used to explore the general impact of cytokines expression on the clustering of the examined COVID-19 patient samples according to disease severity. The result of the unsupervised hierarchical and k-means clustering showed that the collective cytokines panel had little impact on the clustering of the samples according to disease severity, as cases of different degrees of severity were intermingled in both clustering approaches; suggesting an overlap in the signature of some cytokines across the different severity groups. However, the unsupervised clustering gave hints of separation between severe and moderate cases, suggesting that some of the cytokines might have the potential to stratify patients according to disease severity. Therefore, mathematical modeling was carried out to explore further and identify the cytokines that significantly associate with disease severity and other aspects of COVID-19 pathogenesis.

3.4 Optimal Parameter Selection Using Mathematical Modeling

To filter out the biological overlap between the severity groups in the data set and identify key cytokines and biochemical markers that can be used to stratify the clinical parameters collected in the study, an approach combining two mathematical models (multivariate ANOVA with Bonferroni’s stringent multiple testing and Stepwise linear regression) were used (Supplementary Table S2).

3.4.1 Mathematical Modeling Identifies IL-10 as a Biomarker of Severity

Multivariate ANOVA with Bonferroni’s stringent multiple testing was used to determine whether there were statistically significant differences in the expression of particular cytokines and biochemical markers between the COVID-19 severity groups (Figure 4A). The analysis revealed that the levels of IL-10, ANC, ALC, CRP, Ferritin, LDH, BUN, and WBCs were significantly higher in severe cases in comparison to mild-moderate cases. The stepwise linear regression model identified IL-10, PD-L1, TNF -α as potential predictors of COVID-19 disease severity. The data from these two mathematical models suggest a panel of cytokines and biochemical markers for stratifying COVID-19 patients according to disease severity, with the circulating marker IL-10 as the driver key marker.

FIGURE 4

Figure 4 Key driver predictors identified from Multivariate ANOVA with Bonferroni’s stringent multiple testing for (A) disease severity, (B) the requirement for oxygen support, (C) Radiological findings, and (D, E) abnormal liver function indicated by (D) ALT and (E) AST. Means of the predictors’ levels presented as a function of the target variables categories.

3.4.2 Mathematical Modeling Identifies IL-10 as a Biomarker of Oxygen Support Requirement

A similar analysis was performed to determine the potential association between the level of cytokines and biochemical markers and other aspects of COVID-19 pathogenesis, such as the need for oxygen support. Multivariate ANOVA testing suggested the significant association between the need for oxygen support and the levels of GM-CSF, IL-1β, IL-10, ANC, CRP, Ferritin, LDH, BUN, and WBC. Multivariate ANOVA with Bonferroni’s stringent multiple testing revealed the significant increase in the levels of IL-1β and IL-10 in patients requiring mechanical ventilation as opposed to patients depending on room air; and a significant increase in the levels of ANC, LDH, BUN, and WBCs in patients requiring invasive mechanical ventilation in comparison to patients requiring non-invasive forms of oxygen support (e.g., nasal cannula, high flow oxygen mask, and non-invasive positive pressure ventilation) (Figure 4B). The stepwise linear regression model identified IL-10 as a potential predictor of the need for oxygen support. Taken together, both mathematical models suggest IL-10 as a potential marker for the requirement for oxygen support in addition to its potential in stratifying disease severity.

3.4.3 Mathematical Modeling Identifies IL1Ra and IFN-γ as Biomarkers of COVID-19-Specific Radiological Findings

Analysis of the association between chest X-ray (CXR) findings and the levels of cytokines and biochemical markers using multivariate ANOVA with Bonferroni’s stringent multiple testing revealed the upregulated levels of IFN-γ and PD-L1 in normal cases as opposed to patients presenting with consolidation or ground-glass opacities (Figure 4C). On the other hand, IL-1Ra, IL-6, MCP-1, and D-dimer levels were elevated in cases presenting with pneumothorax compared to normal cases or cases presenting with consolidation or ground-glass opacities. Stepwise linear regression analysis proposed IL-1Ra, IFN-γ, and PD-L1 as potential predictors of radiological findings. Stepwise linear regression analysis of CORADs reports suggested IL1Ra and IFN-γ as predictors of radiological findings, further cross-validating the CXR analysis results. Taken together, these data suggest that IL-1Ra and IFN-γ might potentially be used to stratify patients according to radiological findings; IL-1Ra as a potential marker for the development of pneumothorax and IFN-γ as a potential marker predicting the absence of COVID-19 related chest abnormalities.

3.4.4 Mathematical Modeling Identifies IL-6 and Granzyme B as Biomarkers of Liver Injury and Dysfunction

The stepwise linear regression model identified IL-6 and granzyme B as potential predictors of liver injury and dysfunction (indicated by an elevation in the levels of ALT and/or AST). IL-6 and granzyme B levels were elevated in cases with abnormal ALT levels (Figure 4D), while IL-6 was elevated in patients with abnormal levels of AST (Figure 4E).

Intriguingly, Multivariate ANOVA with Bonferroni’s multiple testing revealed the significant reduction in the protein expression level of IL-6 in patients that received COVID-19 treatments (e.g., tocilizumab, lopinavir/ritonavir, favipiravir) in comparison to untreated patients. Treatment status associated significantly with reduced levels of other markers, including IL-1Ra, MCP-1, PD-L1, ALT, D-Dimer, and Albumin.

3.5 Validation of Predictor-Based Stratification of Severity Groups in COVID-19

Cytokines proposed as predictors of disease severity by the two mathematical models were used for model reduction to enhance patients’ clustering. The supervised hierarchical and k-means clustering revealed an enhanced clustering of patients according to disease severity (Figure 5A, B); where severe cases were enriched in the cluster indicated by the red brackets in the heat map (Figure 5A) and clusters 1, 4, and 6 in the k-means PCA plot (Figure 5B). Moreover, ROC curve analysis was used to assess the predictive efficacy of the predictors identified using the mathematical modeling approach to stratify patients according to disease severity. Analysis of the collectively identified cytokines from the two mathematical models (IL-1-α, IL-4, IL-10, IL-13, PD-L1, TNF-α) revealed a significant predictive efficacy with an area under the curve (AUC) value of 0.935. Similarly, assessment of the predictive capacity of the collective biochemical markers identified using multivariate analysis (ANC, Ferritin, LDH, BUN, and WBC) confirmed a significant predictive efficacy with an AUC value of 0.981, Supplementary Figure S2.

FIGURE 5

Figure 5 (A) Heat map representation of the unsupervised hierarchical clustering and (B) Principal Component Analysis (PCA) plot representation of the k-means clustering analysis of cytokines protein expression in the blood samples of COVID-19 patients of different degrees of severity (3 mild, 17 moderate, and 17 severe). ROC analysis of the predictive capacity of the cytokines (AUC=0.93 ± 0.037, 95% CI=0.86-1, p<0.0001). ROC analysis of the predictive capacity of the biochemical markers (AUC=0.98 ± 0.02, 95% CI=0.94-1, p<0.0001), identified using the mathematical models to stratify COVID-19 patients according to disease severity.

ROC analysis of each cytokine and biochemical markers was performed to suggest potential cut-off values with high sensitivity and specificity and significantly high predictive efficacy accordingly. The analysis revealed a cut-off value of 204.5 pg/ml for IL-10, 117.27 pg/ml for PD-L1, 724.0 ng/mL for ferritin, 325.0 U/L for LDH, 10.25×10³/μL for WBC, and 28.27mg/dL for BUN. Cutoff values of the identified predictors for other variables are listed in Supplementary Table S3.

3.6 Mapping of Significantly Differentiated Cytokines on the KEGG Pathways

We further mapped the “predictor” cytokines on several immune-related KEGG pathways, as well as the SARS-CoV2 entry pathway. Of interest, several significantly elevated cytokines in severe COVID-19 patients are remarkable key players along Natural Killer (NK) cell-mediated cytotoxicity pathway (Supplementary Figure S3). Figure 6 summarizes the workflow and the main results.

FIGURE 6

Figure 6 Graphical Abstract of the work flow and the main results.

4 Discussion

In this study, we aimed to predict the outcome of COVID-19 using a non-linear mathematical model of serum cytokine changes, in addition to clinical, biochemical, and radiological parameters. Compared to previous studies, we used AI techniques to integrate data from different modalities (clinical parameters, biochemical tests, cytokine assays, and radiological data) for the first time. We included 87 parameters as input to our model, covering 24 cytokines classified as pro-inflammatory, anti-inflammatory, chemokines, checkpoint markers, receptors and cytotoxic mediators. Cytokines were selected based on an initial transcriptomics analysis of nasopharyngeal swabs of COVID-19 patients and control subjects. Although the unsupervised hierarchical and k-means clustering showed that the collective cytokines panel had little impact on the clustering of the samples, the supervised clustering gave hints of separation between severe and moderate cases, identifying the cytokines with potential predictive value for COVID-19 severity. Taking the initial large number of clinical parameters, biochemical markers and cytokines expression data that would result in a vast number of permutations, machine learning and mathematical modeling were used to filter the data to achieve model reduction and identify the optimum associations to stratify patients according to multiple aspects of COVID-19 pathogenesis. Given the emerging new variants of the virus, our modeling strategy can be applied to different datasets to predict outcomes in new cases of COVID-19 and identify fundamental immune-mediated mechanisms and potential therapeutic targets for such new variants.

Interestingly, in our study, there was a significant upregulation of IL-10 and IL-15, consistently associated with disease severity in both investigated COVID-19 whole blood dataset and our cytokine assays, suggesting the potential utility of these predictive cytokines as circulating biomarkers of severity. IL-10 was reported to contribute to the suppression of the immune system, viral control, and disease severity (31), and a predictor of poor outcomes in COVID-19 patients (32–34). This was possibly linked to its role as an anti-inflammatory cytokine, released as negative feedback in response to the rapid accumulation of pro-inflammatory cytokines (33, 34), and aids in alleviating the CS and preventing tissue damage (31).

Therefore, recombinant IL-10 has been suggested by some investigators for treating acute respiratory distress syndrome (ARDS) in COVID-19 patients based on its immune-regulatory and anti-fibrotic functions (32). Moreover, IL-15/IL-15R axis plays a pivotal role in the function of NK cells (35). IL-15 is produced by activated monocytes/macrophages and activates human NK cells through components of the IL-2R in a pattern similar to that of IL-2. IL-15 also induces IL-10 expression by the NK cells, enhancing its cytotoxic effect. The effect of IL-10 on NK cells is mediated through STAT3 signaling according to in-vitro studies (36). Furthermore, Wang et al., 2021, deciphered that IL-10 regulates metabolic reprogramming in NK cells, via stimulation of the mammalian target of rapamycin complex 1 (mTORC1). In that way, it upregulates glycolysis as well as oxidative phosphorylation in NK cells, thus might enhance the functions of the NK cells (37). Masselli E. et al. reported that the IL-15/IL-15R axis was among the top pathways associated with severe/fatal disease in the viral pandemic gene signature. This may shed some light on the mechanistic role of this axis in the immune NK cell response derangement, which leads to NK cells exhaustion, senescence, apoptosis, and viral persistence (38). Interestingly, severe COVID-19 resulted in an increase of “armed” NK cells containing high levels of cytotoxic proteins such as perforin (39). NK cells are obviously major players in the immune response during COVID-19 infection, but similar to hepatitis virus infections, they may become dysfunctional during severe disease, and their role in organ dysfunction (e.g., liver) requires further investigation (40). It was reported that NK cells might undergo pyroptosis, releasing inflammatory cytokines such as IL-1β (41). However, it is not clear yet whether elevation of IL-15 may activate NK cells, thus contributing to the cytokine burst observed in these patients. In addition, IL-15 was previously suggested to play a role in granulomatous pulmonary diseases through its stimulation of Th-1-driven inflammation (42). Recently, it was also reported to be involved in the development of rapidly progressive interstitial lung disease in polymyositis/dermatomyositis (43). Although reported to enhance NK-cell cytotoxicity, IL-10 elevation may be also a consequence of NK-cell stimulation in an attempt to ameliorate the prevalent hyper-inflammatory state of the CS.

The role of NK cells in COVID-19 infections has not been examined in detail, although it was suggested that they may be important participants (44). In this study we observed that the level of IL-15, a vital cytokine for NK-cell activity, is highly increased, suggesting that NK cells might play a role. It was previously reported that these cells secrete a variety of inflammatory cytokines and chemokines (45), which may contribute to the cytokine storm described in COVID-19 patients.

Several patterns of cytokine changes were identified across the severity levels of COVID-19. Our results showed significant changes in pro-inflammatory cytokines (GM-CSF, IL-6, IL-15, IFN-α, TNF-α, IL-17A) that play a crucial role in the CS (46). The SARS-CoV-2 infection causes local innate immune cells to produce such inflammatory cytokines upon infection of the respiratory epithelial tissue and cause the activation of the adaptive immune cells leading to respiratory epithelial damage (47). This is further supported by activating the inflammasome and NF-κB pathways, inducing the stimulation of several pro-inflammatory genes and immune cell hyperactivation, thus boosting systemic inflammation. In the setting of inflammation, IL-6, which is generated by macrophages and dendritic cells, is known to be a key activator of the JAK/STAT3 pathway (48). Also, IL-6 was reported to contribute to immune cell hyperactivation and target organ dysfunction in COVID-19. On the other hand, the destruction of epithelial cells in the alveolar space caused by SARS-CoV-2 triggers macrophages hyperactivation leading to the CS. IL-6 was found to suppress T lymphocyte activation, that could contribute to lymphopenia in COVID-19 patients. Also, low numbers of T lymphocytes were observed with ICU patients showing high IL-6 and TNF-α serological levels (49–51). Similarly, IL-17 was found to exacerbate lung injury and decrease the survival through the recruitment of neutrophils and stimulation of pro-inflammatory factors (52). GM-CSF also triggers myelopoiesis in order to recruit myeloid cells to the inflammatory sites (53). It was previously suggested as a potential therapy for the COVID-19 CS (47).

Noteworthy, the concurrent elevations in IL-10 and various pro-inflammatory cytokines, and the observed relationship between elevated IL-10 levels and disease severity, suggest that IL-10 is either failing to appropriately suppress inflammation (as observed in other inflammatory conditions (54, 55) or acting in a manner that deviates from its traditional role as an anti-inflammatory molecule, indicating the ability of IL-10 to have different functions under different conditions (32).

The anti-inflammatory IL-4, along with IL-13, mediate the Th2 cell response and M2 polarization, leading to consequent fibrosis and release of growth factors, such as transforming growth factor-β and platelet-derived factor (46, 56). IL-1Ra is known to control the inflammatory immune response by binding to the IL-1R and regulating the production of inflammatory cytokines such as IL-1 and TNF-α (57). In COVID-19 infection, IL-1Ra was suggested to affect the stimulation of pro-inflammatory and antiviral cytokines, where its high level could be an indication of an overactive immune response, thus leading to inflammation-induced tissue damage (34). Controversial patterns were reported regarding IL-4 in previous studies where some reported an increase in peripheral blood/serum of severe COVID-19 patients (30, 46, 58), while others claimed that it did not show any difference (34, 59).

We observed increased levels of chemokines MCP-1 (CCL2), IP-10 (CXCL10), MIP1β (CCL4) and IL-8 (CXCL8) in COVID-19 patients. These chemokines are known to be crucial contributors to pulmonary pathogenesis, such as that observed in COVID-19. CCL2 is known to be released by alveolar macrophages, T cells and endothelial cells in order to induce the migration of inflammatory monocytes and neutrophils along with procollagen synthesis by fibroblasts (60). CXCL10, a known chemoattractant for monocytes, NK and T cells (61), was also reported to play a crucial role in pulmonary neutrophil infiltration (62). Moreover, CXCL10/CXCR3 axis triggers the oxidative burst which promotes exacerbation of the pulmonary inflammation and progression to ARDS (62). CCL4 acts through CCR5 receptor to attract macrophages, dendritic cells, NK and T cells to the site of inflammation (63). Interestingly, the CCR5 antagonist, maraviroc (an antiretroviral medication) was repurposed for moderate to severe COVID-19 (NCT04435522 and NCT04441385). CXCL8 was reported to be responsible for the recruitment, activation, and accumulation of neutrophils (64). Furthermore, it induces the formation of neutrophil extracellular traps (NETs) that further promote inflammation and tissue injury (65). Elevated CXCL8 levels at the time of hospitalization, along with IL-6 and TNF-α, was previously suggested as strong and independent predictors of survival in COVID-19 (20).

PD-L1 was higher in severe COVID-19 patients, whereas the levels of CD40 ligand and granzyme B showed a significant increase in mild-moderate COVID-19 patients, but reduced in severe patients. PD-L1 induces inhibitory signals and apoptosis of CD8+ T cells. This is induced by binding of the pro-inflammatory cytokines to their respective receptors. Hence, the release of IL-6, IL-17, and TNF-α, along with the increased activity of macrophages and neutrophils, cause the increased expression of PD-L1 on the surfaces of immune cells in COVID-19 (66), through the STAT3, PI3K/Akt and NF-κB pathways. CD40L, a costimulatory molecule present on T cells, was found to be released by activated platelets in the serum, that may contribute to pulmonary thrombotic complications observed in COVID-19 as well as being associated with ARDS status (67, 68). Previously, studies have shown that PD-L1 could be a potential predictive factor in various types of cancer (69). The cytotoxic mediator, granzyme B, along with perforin, are the main mediators through which NK cells and cytotoxic T lymphocytes eliminate virally infected host cells, as in COVID-19 infection (70). Interestingly, NK cells from COVID-19 patients exhibit higher levels of granzyme B that is associated with the severity of the disease (71). To re-iterate, our study emphasizes the role of NK cells in COVID-19 infection, an enigma that was not previously resolved. In addition, previous reports revealed that IL-10 and PD-L1 suppress T-cell activity during persistent viral infection (72), thus giving mechanistic insight towards persistent COVID-19 and the potential role of targeting both cytokines to minimize the long-term sequelae of the disease.

Added to combining selective recognized biochemical markers of COVID-19 severity (ANC, CRP, LDH, BUN and ferritin), the triad of elevated serum IL-10, PD-L1 and TNF-α improved the current model accuracy to predict the severity of the disease, through the stepwise linear regression model. The results from these two mathematical models suggest the circulating marker IL-10 as a driving key marker for the stratification of COVID-19 patients according to disease severity. Noteworthy, IL-10 is elevated earlier than IL-6 in COVID-19 patients (32, 34).

Our mathematical model identified IL1-α, IL-4 as negative predictors of severity. As both are involved in adaptive immunity, highlighting its marked derangement in severe COVID-19. In contrast to our results, other reports showed elevated IL1-α in severe COVID-19 that was strongly associated with lung injury (preprint by Liu et al., 2020). Controversial patterns were reported regarding IL-4 plasma levels where some reported an increase in peripheral blood/serum of severe COVID-19 patients (30, 46, 58, 73), while studies show any difference (34, 59).

IL1Ra and IFN-γ were identified in our model, as biomarkers of COVID-19-specific radiological findings. The IL-1 superfamily was previously recognized as a key mediator of inflammation and fibrosis in different organs, with IL-1Ra as an antagonistic cytokine (74). The crucial balance between IL-1β and IL1Ra determines the resultant immune response in many tissues (74). In the severe COVID-19 cases in this study, IL-1β significantly decreased and IL1Ra significantly increased as a part of the marked immune dysregulation. This was associated with specific COVID-19 related radiological findings, as revealed by the mathematical model. IFN-γ mediates immune-mediated damage in acute lung injury (75).

In support of our findings related to IL-6 and granzyme B as biomarkers of liver injury, a recent study demonstrated that IL-6 trans-signaling drives COVID-19-associated hepatic endotheliopathy, which is suggested as a possible mechanism underlying the liver injury (29). Previous reports highlighted the role of NK cells and their enzymes (Granzyme B and perforin) in hepatic immune homeostasis (76). IL-6 was reported to suppress the NK cytotoxicity in-vitro and in-vivo (77). However, in view of the multiple cytokines affecting the NK cells in the CS context and the elevation of IL-15 (NK stimulator), the effect of IL-6 is surpassed, with a net result of increased granzyme B. Our model shows the high accuracy of liver injury prediction in severe COVID-19, by combining IL-6 and granzyme B as predictors.

We used a stochastic non-linear modeling approach to reduce the dataset for multi-dimensional data and to integrate data from different modalities. The non-linear ODE model is crucial to clearly reflect the dynamics of biological systems (78). To estimate the exact probabilities for biological systems, approaches are mainly based on Monte Carlo sampling (e.g. the Stochastic Simulation Algorithm) (79). To create a dynamic model of CS, Waito et al., 2016 used a nonlinear differential equation model, considering the cytokine production rate in relation to their interactions with one another. They adjusted the model by using the data from a CS mouse model (IFN type 1 receptor KO). Interestingly and concordant to our results, the model revealed that TNF-α, IL-10, IL-6, and MIP-1β, exerted the largest effects on the dynamics of the cytokine storm (17). In the current study, we used non-linear modeling that attempts to identify global solutions to integrate and explore biomarkers that can predict COVID-19 severity (8).

Our study sheds light on key immunological aspects of the COVID-19-CS that seem to significantly differ from the CS occurring in other diseases. Beyond its value as a biological predictive tool, our mathematical analysis poses important questions for future research.

5 Conclusions

Predictive modeling in COVID-19 has gained a high value, considering the complexity of the disease. Using a non-linear model for clinical, biochemical, immunological, and radiological data could achieve a high level of prediction accuracy. In our proposed integrative model, we validated a cytokine panel derived from transcriptomics analysis of nasopharyngeal swab samples of COVID-19 patients. Our model advocates the trio of IL-10, PD-L1 and TNF-α as an accurate predictor of severity, in addition to previously recognized ANC, CRP, LDH, BUN and ferritin, whereas IL-1α, IL-4 were negative predictors. IL-10 was shown to be a driving marker and a positive predictor of mechanical ventilation. Moreover, IFN-γ, IL-1Ra were predictors of remarkable radiological findings, whereas high IL-6 and granzyme B were found to predict liver injury in COVID-19 patients.

We identified key cytokines that were consistently associated with severity, like IL-10, an enhancer of NK cytotoxicity, and IL-15, a stimulator of NK cells, Obviously, the modeling methodology can be used to identify key players and predict outcome in new variants of COVID-19.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://figshare.com/articles/dataset/Transcriptome_COVID19_Cytokine_Mathmodel_xlsx/19386194.

Ethics Statement

The studies involving human participants were reviewed and approved by Approval of the Research and Ethics Committee at Dubai Health Authority (DSREC-04/2020_09) and (DSREC-04/2020_19). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

Conceptualization: MS-A and RH. Data curation: NE, SH, LS, and MS. Formal analysis: NE, SH, PH, and MS. Funding acquisition: MS-A and RH. Investigation: NE, SH, and IT. Methodology: MS-A, RH, NE, and SH. Project administration: MS-A. Resources: LS, BM, HS, MS, KOH, OA-A, JT, NS, AM, and QH. Software and mathematical modeling: SH. Supervision: RH, BM, HS, and MS-A, Validation: NE. Visualization: NE, SH, and PH. Writing - original draft: NE, SH, IT, MS-A, JT, and RH. Writing - review and editing: HS, AM, and QH. All authors contributed to the article and approved the submitted version.

Funding

MS-A and RH are funded by the University of Sharjah Research Grants CoV-19 #0304 and CoV-19 #0308.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to acknowledge the efforts and help of the COVID-lab team at the University of Sharjah, Prof. Rabih Halwani, Dr Abdul Wahid Ansari, Dr Narjes Saheb Sharif-Askari, and Dr Fatemeh Saheb Sharif-Askari.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.865845/full#supplementary-material

References

1. Fajgenbaum DC, June CH. Cytokine Storm. N Engl J Med (2020) 383(23):2255–73. doi: 10.1056/NEJMra2026131

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Xu Z-S, Shu T, Kang L, Wu D, Zhou X, Liao B-W, et al. Temporal Profiling of Plasma Cytokines, Chemokines and Growth Factors From Mild, Severe and Fatal COVID-19 Patients. Signal Transduct Targeted Ther (2020) 5(1):100. doi: 10.1038/s41392-020-0211-1

CrossRef Full Text | Google Scholar

3. Guo J, Wang S, Xia H, Shi D, Chen Y, Zheng S, et al. Cytokine Signature Associated With Disease Severity in COVID-19. Front Immunol (2021) 12:681516. doi: 10.3389/fimmu.2021.681516

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ramatillah DL, Gan SH, Pratiwy I, Syed Sulaiman SA, Jaber AAS, Jusnita N, et al. Impact of Cytokine Storm on Severity of COVID-19 Disease in a Private Hospital in West Jakarta Prior to Vaccination. PloS One (2022) 17(1):e0262438-e. doi: 10.1371/journal.pone.0262438

CrossRef Full Text | Google Scholar

5. Furlow B. COVACTA Trial Raises Questions About Tocilizumab’s Benefit in COVID-19. Lancet Rheumatol (2020) 2(10):e592. doi: 10.1016/S2665-9913(20)30313-1

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Telenti A, Arvin A, Corey L, Corti D, Diamond MS, García-Sastre A, et al. After the Pandemic: Perspectives on the Future Trajectory of COVID-19. Nature (2021) 596(7873):495–504. doi: 10.1038/s41586-021-03792-w

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ioannidis J, Cripps S, Tanner MA. Forecasting for COVID-19 has Failed. Int Institute Forecasters (2020) 38(2):423–38. doi: 10.1016/j.ijforecast.2020.08.004

CrossRef Full Text | Google Scholar

8. Raue A, Schilling M, Bachmann J, Matteson A, Schelke M, Kaschek D, et al. Lessons Learned From Quantitative Dynamical Modeling in Systems Biology. PloS One (2013) 8(9):e74335. doi: 10.1371/journal.pone.0074335

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Xavier MP, Bonin CRB, dos Santos RW, Lobosco M eds. International Conference on Bioinformatics and Biomedicine (BIBM). Kansas City, MO: IEEE (2017).

Google Scholar

10. Almocera AES, Quiroz G, Hernandez-Vargas EA. Stability Analysis in COVID-19 Within-Host Model With Immune Response. Commun Nonlinear Sci Numer Simul (2021) 95:105584. doi: 10.1016/j.cnsns.2020.105584

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Du SQ, Yuan W. Mathematical Modeling of Interaction Between Innate and Adaptive Immune Responses in COVID-19 and Implications for Viral Pathogenesis. J Med Virol (2020) 92(9):1615–28. doi: 10.1002/jmv.25866

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Hernandez-Vargas EA, Velasco-Hernandez JX. In-Host Mathematical Modelling of COVID-19 in Humans. Annu Rev Control (2020) 50:448–56. doi: 10.1016/j.arcontrol.2020.09.006

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xavier MP, Reis RF, dos Santos RW, Lobosco M eds. International Conference on Bioinformatics and Biomedicine (BIBM). Seoul, South Korea: IEEE (2020).

Google Scholar

14. Mahboub B, Bataineh MTA, Alshraideh H, Hamoudi R, Salameh L, Shamayleh A. Prediction of COVID-19 Hospital Length of Stay and Risk of Death Using Artificial Intelligence-Based Modeling. Front Med (2021) 8:389. doi: 10.3389/fmed.2021.592336

CrossRef Full Text | Google Scholar

15. Shiri I, Sorouri M, Geramifar P, Nazari M, Abdollahi M, Salimi Y, et al. Machine Learning-Based Prognostic Modeling Using Clinical Data and Quantitative Radiomic Features From Chest CT Images in COVID-19 Patients. Comput Biol Med (2021) 132:104304. doi: 10.1016/j.compbiomed.2021.104304

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Wang JM, Liu W, Chen X, McRae MP, McDevitt JT, Fenyö D. Predictive Modeling of Morbidity and Mortality in Patients Hospitalized With COVID-19 and Its Clinical Implications: Algorithm Development and Interpretation. J Med Internet Res (2021) 23(7):e29514. doi: 10.2196/29514

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Waito M, Walsh SR, Rasiuk A, Bridle BW, Willms AR eds. A Mathematical Model of Cytokine Dynamics During a Cytokine Storm. Mathematical and Computational Approaches in Advancing Modern Science and Engineering. Cham: Springer International Publishing (2016).

Google Scholar

18. Yiu HH, Graham AL, Stengel RF. Dynamics of a Cytokine Storm. PloS One (2012) 7(10):e45027. doi: 10.1371/journal.pone.0045027

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Rooney C, Sauer T. Modeling Cytokine Release Syndrome. Nat Med (2018) 24(6):705–6. doi: 10.1038/s41591-018-0068-9

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Del Valle DM, Kim-Schulze S, Huang H-H, Beckmann ND, Nirenberg S, Wang B, et al. An Inflammatory Cytokine Signature Predicts COVID-19 Severity and Survival. Nat Med (2020) 26(10):1636–43. doi: 10.1038/s41591-020-1051-9

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Perreau M, Suffiotti M, Marques-Vidal P, Wiedemann A, Levy Y, Laouénan C, et al. The Cytokines HGF and CXCL13 Predict the Severity and the Mortality in COVID-19 Patients. Nat Commun (2021) 12(1):4888. doi: 10.1038/s41467-021-25191-5

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Penha D, Pinto EG, Matos F, Hochhegger B, Monaghan C, Taborda-Barata L, et al. CO-RADS: Coronavirus Classification Review. J Clin Imaging Sci (2021) 11:9. doi: 10.25259/JCIS_192_2020

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Hammoudeh SM, Venkatachalam T, Ansari AW, Bendardaf R, Hamid Q, Rahmani M, et al. Systems Immunology Analysis Reveals an Immunomodulatory Effect of Snail-P53 Binding on Neutrophil- and T Cell-Mediated Immunity in KRAS Mutant Non-Small Cell Lung Cancer. Front Immunol (2020) 11:569671. doi: 10.3389/fimmu.2020.569671

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc (Methodological) (1995) 57(1):289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

25. Bernardes JP, Mishra N, Tran F, Bahmer T, Best L, Blase JI, et al. Longitudinal Multi-Omics Analyses Identify Responses of Megakaryocytes, Erythroid Cells, and Plasmablasts as Hallmarks of Severe COVID-19. Immunity (2020) 53(6):1296–314.e9. doi: 10.1016/j.immuni.2020.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wei C, Li J, Bumgarner RE. Sample Size for Detecting Differentially Expressed Genes in Microarray Experiments. BMC Genomics (2004) 5:87. doi: 10.1186/1471-2164-5-87

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Trevethan R. Sensitivity, Specificity, and Predictive Values: Foundations, Pliabilities, and Pitfalls in Research and Practice. Front Public Health (2017) 5:307. doi: 10.3389/fpubh.2017.00307

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ziegler CGK, Miao VN, Owings AH, Navia AW, Tang Y, Bromley JD, et al. Impaired Local Intrinsic Immunity to SARS-CoV-2 Infection in Severe COVID-19. Cell (2021) 184(18):4713–33.e22. doi: 10.1016/j.cell.2021.07.023

PubMed Abstract | CrossRef Full Text | Google Scholar

29. McConnell MJ, Kawaguchi N, Kondo R, Sonzogni A, Licini L, Valle C, et al. Liver Injury in COVID-19 and IL-6 Trans-Signaling-Induced Endotheliopathy. J Hepatol (2021) 75(3):647–58. doi: 10.1016/j.jhep.2021.04.050

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, et al. Clinical Features of Patients Infected With 2019 Novel Coronavirus in Wuhan, China. Lancet (2020) 395(10223):497–506. doi: 10.1016/S0140-6736(20)30183-5

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Islam H, Chamberlain TC, Mui AL, Little JP. Elevated Interleukin-10 Levels in COVID-19: Potentiation of Pro-Inflammatory Responses or Impaired Anti-Inflammatory Action? Front Immunol (2021) 12:2485. doi: 10.3389/fimmu.2021.677008

CrossRef Full Text | Google Scholar

32. Lu L, Zhang H, Dauphars DJ, He YW. A Potential Role of Interleukin 10 in COVID-19 Pathogenesis. Trends Immunol (2021) 42(1):3–5. doi: 10.1016/j.it.2020.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Han H, Ma Q, Li C, Liu R, Zhao L, Wang W, et al. Profiling Serum Cytokines in COVID-19 Patients Reveals IL-6 and IL-10 Are Disease Severity Predictors. Emerg Microbes Infect (2020) 9(1):1123–30. doi: 10.1080/22221751.2020.1770129

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhao Y, Qin L, Zhang P, Li K, Liang L, Sun J, et al. Longitudinal COVID-19 Profiling Associates IL-1RA and IL-10 With Disease Severity and RANTES With Mild Disease. JCI Insight (2020) 5(13):e139834. doi: 10.1172/jci.insight.139834

CrossRef Full Text | Google Scholar

35. Fehniger TA, Suzuki K, VanDeusen JB, Cooper MA, Freud AG, Caligiuri MA. Fatal Leukemia in Interleukin-15 Transgenic Mice. Blood Cells Mol Dis (2001) 27(1):223–30. doi: 10.1006/bcmd.2001.0379

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Park JY, Lee SH, Yoon S-R, Park Y-J, Jung H, Kim T-D, et al. IL-15-Induced IL-10 Increases the Cytolytic Activity of Human Natural Killer Cells. Mol Cells (2011) 32(3):265–72. doi: 10.1007/s10059-011-1057-8

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Wang Z, Guan D, Huo J, Biswas SK, Huang Y, Yang Y, et al. IL-10 Enhances Human Natural Killer Cell Effector Functions via Metabolic Reprogramming Regulated by Mtorc1 Signaling. Front Immunol (2021) 12:619195. doi: 10.3389/fimmu.2021.619195

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Masselli E, Vitale M. NK Cells on the ViP Stage of COVID-19. EBioMedicine (2021) 69:103458. doi: 10.1016/j.ebiom.2021.103458

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Maucourant C, Filipovic I, Ponzetta A, Aleman S, Cornillet M, Hertwig L, et al. Natural Killer Cell Immunotypes Related to COVID-19 Disease Severity. Sci Immunol (2020) 5(50):eabd6832. doi: 10.1126/sciimmunol.abd6832

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Highton AJ, Schuster IS, Degli-Esposti MA, Altfeld M. The Role of Natural Killer Cells in Liver Inflammation. Semin Immunopathol (2021) 43(4):519–33. doi: 10.1007/s00281-021-00877-6

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hachim MY, Khalil BA, Elemam NM, Maghazachi AA. Pyroptosis: The Missing Puzzle Among Innate and Adaptive Immunity Crosstalk. J Leukocyte Biol (2020) 108(1):323–38. doi: 10.1002/JLB.3MIR0120-625R

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Hamid Q, Ito I, Muro S. Laurent GJ, Shapiro SD, editors. Encyclopedia of Respiratory Medicine. Oxford: Academic Press (2006). p. 385–90.

Google Scholar

43. Schnabel A, Reuter M, Biederer J, Richter C, Gross WL. Interstitial Lung Disease in Polymyositis and Dermatomyositis: Clinical Course and Response to Treatment. Semin Arthritis Rheum (2003) 32(5):273–84. doi: 10.1053/sarh.2002.50012

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Al-Ani M, Elemam NM, Hundt JE, Maghazachi AA. Drugs for Multiple Sclerosis Activate Natural Killer Cells: Do They Protect Against COVID-19 Infection? Infect Drug Resist (2020) 13:3243–54. doi: 10.2147/IDR.S269797

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Rolin J, Sand KL, Knudsen E, Maghazachi AA. FTY720 and SEW2871 Reverse the Inhibitory Effect of S1P on Natural Killer Cell Mediated Lysis of K562 Tumor Cells and Dendritic Cells But Not on Cytokine Release. Cancer Immunol Immunother CII (2010) 59(4):575–86. doi: 10.1007/s00262-009-0775-7

CrossRef Full Text | Google Scholar

46. Mehta P, McAuley DF, Brown M, Sanchez E, Tattersall RS, Manson JJ, et al. COVID-19: Consider Cytokine Storm Syndromes and Immunosuppression. Lancet (2020) 395(10229):1033–4. doi: 10.1016/S0140-6736(20)30628-0

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Yang L, Xie X, Tu Z, Fu J, Xu D, Zhou Y. The Signal Pathways and Treatment of Cytokine Storm in COVID-19. Signal Transduct Targeted Ther (2021) 6(1):255. doi: 10.1038/s41392-021-00679-0

CrossRef Full Text | Google Scholar

48. Kang S, Tanaka T, Narazaki M, Kishimoto T. Targeting Interleukin-6 Signaling in Clinic. Immunity (2019) 50(4):1007–23. doi: 10.1016/j.immuni.2019.03.026

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Sarzi-Puttini P, Giorgi V, Sirotti S, Marotto D, Ardizzone S, Rizzardini G, et al. COVID-19, Cytokines and Immunosuppression: What can We Learn From Severe Acute Respiratory Syndrome? Clin Exp Rheumatol (2020) 38(2):337–42.

PubMed Abstract | Google Scholar

50. Ye Q, Wang B, Mao J. The Pathogenesis and Treatment of the `Cytokine Storm’ in COVID-19. J Infect (2020) 80(6):607–13. doi: 10.1016/j.jinf.2020.03.037

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Diao B, Wang C, Tan Y, Chen X, Liu Y, Ning L, et al. Reduction and Functional Exhaustion of T Cells in Patients With Coronavirus Disease 2019 (COVID-19). Front Immunol (2020) 11:827. doi: 10.3389/fimmu.2020.00827

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Muir R, Osbourn M, Dubois AV, Doran E, Small DM, Monahan A, et al. Innate Lymphoid Cells Are the Predominant Source of IL-17A During the Early Pathogenesis of Acute Respiratory Distress Syndrome. Am J Respir Crit Care Med (2016) 193(4):407–16. doi: 10.1164/rccm.201410-1782OC

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Hamilton JA. GM-CSF in Inflammation. J Exp Med (2020) 217(1):e20190945. doi: 10.1084/jem.20190945

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Antoniv TT, Ivashkiv LB. Dysregulation of Interleukin-10–Dependent Gene Expression in Rheumatoid Arthritis Synovial Macrophages. Arthritis Rheumatism (2006) 54(9):2711–21. doi: 10.1002/art.22055

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Neidhart M, Jüngel A, Ospelt C, Michel BA, Gay RE, Gay S. Deficient Expression of Interleukin-10 Receptor α Chain in Rheumatoid Arthritis Synovium: Limitation of Animal Models of Inflammation. Arthritis Rheumatism (2005) 52(10):3315–8. doi: 10.1002/art.21274

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Thompson BT, Chambers RC, Liu KD. Acute Respiratory Distress Syndrome. N Engl J Med (2017) 377(6):562–72. doi: 10.1056/NEJMra1608077

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Oleksowicz L, Dutcher JP. A Review of the New Cytokines: IL-4, IL-6, IL-11, and IL-12. Am J Ther (1994) 1(2):107–15. doi: 10.1097/00045391-199408000-00002

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Lucas C, Wong P, Klein J, Castro TBR, Silva J, Sundaram M, et al. Longitudinal Analyses Reveal Immunological Misfiring in Severe COVID-19. Nature (2020) 584(7821):463–9. doi: 10.1038/s41586-020-2588-y

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Kang CK, Han G-C, Kim M, Kim G, Shin HM, Song K-H, et al. Aberrant Hyperactivation of Cytotoxic T-Cell as a Potential Determinant of COVID-19 Severity. Int J Infect Dis (2020) 97:313–21. doi: 10.1016/j.ijid.2020.05.106

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Khalil BA, Elemam NM, Maghazachi AA. Chemokines and Chemokine Receptors During COVID-19 Infection. Comput Struct Biotechnol J (2021) 19:976–88. doi: 10.1016/j.csbj.2021.01.034

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Liu M, Guo S, Stiles JK. The Emerging Role of CXCL10 in Cancer (Review). Oncol Lett (2011) 2(4):583–9. doi: 10.3892/ol.2011.300

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Ichikawa A, Kuba K, Morita M, Chida S, Tezuka H, Hara H, et al. CXCL10-CXCR3 Enhances the Development of Neutrophil-Mediated Fulminant Lung Injury of Viral and Nonviral Origin. Am J Respir Crit Care Med (2013) 187(1):65–77. doi: 10.1164/rccm.201203-0508OC

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Patterson BK, Guevara-Coto J, Yogendra R, Francisco EB, Long E, Pise A, et al. Immune-Based Prediction of COVID-19 Severity and Chronicity Decoded Using Machine Learning. Front Immunol (2021) 12:2520. doi: 10.3389/fimmu.2021.700782

CrossRef Full Text | Google Scholar

64. Peveri P, Walz A, Dewald B, Baggiolini M. A Novel Neutrophil-Activating Factor Produced by Human Mononuclear Phagocytes. J Exp Med (1988) 167(5):1547–59. doi: 10.1084/jem.167.5.1547

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Saffarzadeh M, Juenemann C, Queisser MA, Lochnit G, Barreto G, Galuska SP, et al. Neutrophil Extracellular Traps Directly Induce Epithelial and Endothelial Cell Death: A Predominant Role of Histones. PloS One (2012) 7(2):e32366. doi: 10.1371/journal.pone.0032366

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Chen J, Vitetta L. Increased PD-L1 Expression May Be Associated With the Cytokine Storm and CD8+ T-Cell Exhaustion in Severe COVID-19. J Infect Dis (2021) 223(9):1659–60. doi: 10.1093/infdis/jiab061

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Poissy J, Goutay J, Caplan M, Parmentier E, Duburcq T, Lassalle F, et al. Pulmonary Embolism in Patients With COVID-19. Circulation (2020) 142(2):184–6. doi: 10.1161/CIRCULATIONAHA.120.047430

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Petrey AC, Qeadan F, Middleton EA, Pinchuk IV, Campbell RA, Beswick EJ. Cytokine Release Syndrome in COVID-19: Innate Immune, Vascular, and Platelet Pathogenic Factors Differ in Severity of Disease and Sex. J Leukocyte Biol (2021) 109(1):55–66. doi: 10.1002/JLB.3COVA0820-410RRR

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Yang Y, Pang Z, Ding N, Dong W, Ma W, Li Y, et al. The Efficacy and Potential Predictive Factors of PD-1/PD-L1 Blockades in Epithelial Carcinoma Patients: A Systematic Review and Meta Analysis. Oncotarget (2016) 7(45):74350–61. doi: 10.18632/oncotarget.11291

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Cunningham L, Kimber I, Basketter D, Simmonds P, McSweeney S, Tziotzios C, et al. Perforin, COVID-19 and a Possible Pathogenic Auto-Inflammatory Feedback Loop. Scandinavian J Immunol (2021) 94(5):e13102. doi: 10.1111/sji.13102

CrossRef Full Text | Google Scholar

71. Zenarruzabeitia O, Astarloa-Pando G, Terrén I, Orrantia A, Pérez-Garay R, Seijas-Betolaza I, et al. T Cell Activation, Highly Armed Cytotoxic Cells and a Shift in Monocytes CD300 Receptors Expression Is Characteristic of Patients With Severe COVID-19. Front Immunol (2021) 12:620. doi: 10.3389/fimmu.2021.655934

CrossRef Full Text | Google Scholar

72. Brooks DG, Ha S-J, Elsaesser H, Sharpe AH, Freeman GJ, Oldstone MBA. IL-10 and PD-L1 Operate Through Distinct Pathways to Suppress T-Cell Activity During Persistent Viral Infection. Proc Natl Acad Sci (2008) 105(51):20428. doi: 10.1073/pnas.0811139106

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Chen N, Zhou M, Dong X, Qu J, Gong F, Han Y, et al. Epidemiological and Clinical Characteristics of 99 Cases of 2019 Novel Coronavirus Pneumonia in Wuhan, China: A Descriptive Study. Lancet (2020) 395(10223):507–13. doi: 10.1016/S0140-6736(20)30211-7

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Borthwick LA. The IL-1 Cytokine Family and Its Role in Inflammation and Fibrosis in the Lung. Semin Immunopathol (2016) 38(4):517–34. doi: 10.1007/s00281-016-0559-z

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Mock JR, Tune MK, Dial CF, Torres-Castillo J, Hagan RS, Doerschuk CM. Effects of IFN-γ on Immune Cell Kinetics During the Resolution of Acute Lung Injury. Physiol Rep (2020) 8(3):e14368. doi: 10.14814/phy2.14368

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Mikulak J, Bruni E, Oriolo F, Di Vito C, Mavilio D. Hepatic Natural Killer Cells: Organ-Specific Sentinels of Liver Immune Homeostasis and Physiopathology. Front Immunol (2019) 10:946. doi: 10.3389/fimmu.2019.00946

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Cifaldi L, Prencipe G, Caiello I, Bracaglia C, Locatelli F, De Benedetti F, et al. Inhibition of Natural Killer Cell Cytotoxicity by Interleukin-6: Implications for the Pathogenesis of Macrophage Activation Syndrome. Arthritis Rheumatol (Hoboken NJ) (2015) 67(11):3037–46. doi: 10.1002/art.39295

CrossRef Full Text | Google Scholar

78. Yates A, Chan CC, Callard RE, George AJ, Stark J. An Approach to Modelling in Immunology. Briefings Bioinf (2001) 2(3):245–57. doi: 10.1093/bib/2.3.245

CrossRef Full Text | Google Scholar

79. Tanevski J, Todorovski L, Džeroski S. Learning Stochastic Process-Based Models of Dynamical Systems From Knowledge and Data. BMC Syst Biol (2016) 10(1):30. doi: 10.1186/s12918-016-0273-4

PubMed Abstract | CrossRef Full Text | Google Scholar