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/m2) 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).
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.