Model-interpreted outcomes of artificial neural networks classifying immune biomarkers associated with severe infections in ICU

Introduction Millions of deaths worldwide are a result of sepsis (viral and bacterial) and septic shock syndromes which originate from microbial infections and cause a dysregulated host immune response. These diseases share both clinical and immunological patterns that involve a plethora of biomarkers that can be quantified and used to explain the severity level of the disease. Therefore, we hypothesize that the severity of sepsis and septic shock in patients is a function of the concentration of biomarkers of patients. Methods In our work, we quantified data from 30 biomarkers with direct immune function. We used distinct Feature Selection algorithms to isolate biomarkers to be fed into machine learning algorithms, whose mapping of the decision process would allow us to propose an early diagnostic tool. Results We isolated two biomarkers, i.e., Programmed Death Ligand-1 and Myeloperoxidase, that were flagged by the interpretation of an Artificial Neural Network. The upregulation of both biomarkers was indicated as contributing to increase the severity level in sepsis (viral and bacterial induced) and septic shock patients. Discussion In conclusion, we built a function considering biomarker concentrations to explain severity among sepsis, sepsis COVID, and septic shock patients. The rules of this function include biomarkers with known medical, biological, and immunological activity, favoring the development of an early diagnosis system based in knowledge extracted from artificial intelligence.


Introduction
Sepsis and septic shock are life-threatening syndromes that are associated with dysregulation in the host immune responses to infection (1). They can lead to organ failure and consequently death (2). As an example, in 2017, more than 11 million deaths associated with sepsis were reported worldwide (3), which represented a mortality rate of approximately 22%. Moreover, a considerable share of viral sepsis patients (i.e., sepsis COVID) meet the definition for sepsis-3 (bacteria-induced sepsis) (4). A more severe manifestation of sepsis is septic shock, in which patients meet all sepsis-3 criteria and require the use of vasopressor (1).
The fact that sepsis (of viral and bacterial sources) and its severe subset (i.e., septic shock) meet the criteria for sepsis-3 definition allows the employment of biomarkers as an early diagnostic tool (5). We selected a consortium composed of the following 30 biomarkers that are responsible for reflecting signals of specific moments during an immune response towards a pathogen that causes sepsis and/or septic shock: The inflammatory pathway of the diseases explored here leaves traces behind that might be employed in profiling the severity of the diseases themselves (6,7). Many of these traces are depicted through the analysis of biomarkers (i.e., cytokines and chemokines) associated with the host immune system. For example, VCAM-1, ICAM-1, and VEGFC are recruited when there is damage to vascular tissue (8). Moreover, CCL2 orchestrates the recruitment of immune cells to sites of inflammation (9). In addition, PDL1 functions as a suppressor of the adaptive immune system as it binds to the Programmed Cell Death Protein 1 (PD1) (10). Finally, MPO is mainly expressed in neutrophil granulocytes, granting antipathogenic activity to the immune cells expressing them (11). Each biomarker has an individual and grouped function in inflammatory pathways (12); therefore, a systematic analysis using data mining to determine the key items involved in sepsis/ septic shock syndrome is appreciated (13,14).
It has been stated (15) that a dataset comprised of too many dimensions (i.e., biomarkers) might slow, mask, and reduce the efficiency of machine learning approaches. Therefore, the selection of the most important features of a dataset is an essential step within data pre-processing frameworks. Moreover, Explainable Artificial Intelligence (XAI) has arisen as a manner to promote comprehension for the decision pattern employed by machine learning approaches, especially with the solid ethical standards required by the medical sciences. Therefore, decision-making processes benefit of a mathematically evidenced procedure (16).
In this paper, we hypothesize that the severity of sepsis (bacterial and viral) and septic shock patients is a function of the concentration of biomarkers. For that, we aim to use feature selection algorithms that will isolate biomarkers as candidates for distinguishing the severity of multi-organ failure in sepsis, sepsis COVID, and septic shock patients. With subsets of the selected biomarkers, we aim to evaluate these subsets through interpretable Artificial Neural Networks, so the biomarker concentration that defines the severity of patients with sepsis, sepsis COVID, and septic shock can be used as an early diagnostic tool.

Study design
All the samples used in this study were obtained from a critically ill cohort of Intensive Care Unit (ICU) sepsis, sepsis COVID, and septic shock patients at St James's Hospital in Dublin, Ireland. Institutional Research Board approval was granted by the SJH/TUH Joint Research Ethics Committee and The Health Research Consent Declaration Committee (HRCDC) under the register number REC: 2020-05 List 17 and project ID 0428. Biological samples, clinical findings, and laboratory data were collected at days 0, 3, and 14 after presentation of severe infection to monitor the progression and sepsis-induced immune-paralysis state at different stages of the disease. Sample collection took place from September 2020 to March 2021. Sequential Organ Failure Assessment (SOFA) score was obtained on admission to the ICU and at the matching collection timepoints for samples. The clinical variables for white blood cell (WBC) count (worst record of day 0), neutrophils (day 0), positive culture, and up to five comorbidities were attributed to each patient.
The plasma concentration of both sIL6R and sGP130 was evaluated and quantified with Enzyme-Linked Immunosorbent Assay (ELISA) kits (BMS214TEN for sIL-6R and EHIL6STX10 for sGP130; ThermoFisher Scientific). These biomarkers were selected due to their crucial role in the IL-6 inflammatory pathway (17). Finally, each sample was assayed and quantified following the kit manufacturer's instructions. All samples were obtained from patients already admitted to the intensive care unit (ICU). We included the concentrations for all biomarkers at day zero of ICU in Supplementary Material S2. We provide the dataset in Supplementary Material S2.
2.3 Binarization of SOFA score into different degrees of multi-organ failure Each patient in our cohort was assigned a Sequential Organ Failure Assessment (SOFA) score. For classification purposes, we binarized the SOFA score into two groups, i.e., High Degree Multi-Organ Failure (HDMOF) and Low Degree Multi-Organ Failure (LDMOF). We employed a cut-off value of 8, as reported by (18) to binarize the groups; thus the HDMOF group is characterized by a SOFA score equal to or higher than 8 while the LDMOF group has a SOFA score less than 8. In their paper, upon binarizing patients with a SOFA score cut of (>=8, and<8) Martin-Loeches et al. (2017) (18) found different mortality rates and antibodies levels that well explained the severity of sepsis patients.

Statistical analyses
All statistical procedures were performed using R. We used the Shapiro-Wilk test to find data distribution. To compare averages between groups, the Kruskal-Wallis test, t test, and ANOVA test were used through the Rstatix package (version 0.7.0). The Kaplan-Meier method under the survival R package (version 3.3-1) was used for calculating survival probabilities for each group (i.e., sepsis, sepsis COVID, and septic shock). The level of significance was set to 0.05.

Feature selection
To systematically choose biomarkers that are associated with high and low MOF, a Feature Selection (FS) step was applied, and the outcomes were benchmarked. The FS algorithms chosen for this study are representatives of three distinct classes of these methods. For a wrapper algorithm, we selected Boruta and applied the parameters specified in (19); the filter algorithm we chose is Information Gain (IG); and a combination of both FS classes, which results in an embedded algorithm has as its representative the Lasso Regression (LR) under the parameters specified in (20). Both Boruta and LR were implemented in R through the packages Boruta (version 7.0.0) and Glmnet (version 4.1-4), respectively. IG was implemented in Python using the class mutual_info_classif found in the sklearn.feature_selection library (version 1.1.0). Both the R and Python scripts that performed the feature selection process are available at https://github.com/gustavsganzerla/covidbiomarker/blob/main/XAI/feature_selection.R and https:// github.com/gustavsganzerla/covid-biomarker/blob/main/XAI/ information_gain.py, respectively.

Classification of severity with artificial intelligence
We selected the algorithms Support Vector Machines (SVM), Random Forest (RF), Classification and Regression Trees (CART), K-Nearest Neighbors (KNN), and deep learning Artificial Neural Networks (ANN) to classify patients' severity with immune biomarkers as input. For that, we performed a 10-fold crossvalidation process on the input data. The test dataset encompassed 20% of the whole data. The first five classification algorithms were implemented in R through the Caret package (version 6.0-9) and their code is available at https://github.com/gustavsganzerla/covidbiomarker/blob/main/XAI/classific.R.
The ANN approach was developed in Python using the Tensorflow library (version 2.8.0). A sigmoid function was used in the output layer to return a probability of a given patient having MOF or not. The outcome was binarized through a confusion matrix built under the default 0.5 decision threshold applied in the outcome of the output neuron. The classes obtained were categorized as follows: True Positives (TPs), True Negatives (TNs), False Positive (FPs), and False Negatives (FNs). The ANN performance was evaluated in terms of accuracy, AUC, FPs, FNs, TPs, TNs, all of which can be found in the module tf.keras.metrics. The architecture of the ANN has input layers that vary according to the set of biomarkers entered. Next, the model has two fully connected layers with 10 hidden neurons each. We also increased the epochs of the ANN until the error (binary cross entropy) kept dropping. The scripts containing the ANN simulation are available at https://github.com/gustavsganzerla/covid-biomarker/blob/main/ XAI/SHAP-ANN.

Explanation of the classification model
We used Shapley Additive Explanations (SHAP) (21) to provide interpretability to the successful classification models. SHAP will assign a score (either positive or negative) for each input variable in assigning a label to an observation. If the SHAP score of a given feature is positive, it is positively correlated with the assigning of a label; otherwise, it is negatively correlated with the target label.
The SHAP approach is defined as a solid theoretical foundation that may be used to explain any predictive model locally and globally. From this, we employed the kernel.explainer method in the SHAP module. As our classifier is not a tree-based algorithm, the Kernel SHAP is applied. Kernel SHAP will measure the contribution of each input feature to the outcome of the model, it consists of five steps: I. Sample all possible coalitions (i.e., combination of input features) in the dataset: (0 = feature absent and 1 = feature present in the coalition). II. The prediction of each Z' K is obtained with the application of Z' K to the predictive model. III. The weight of each Z' K is computed by the SHAP kernel. IV. The model is fitted. V. Return the SHAP coefficients f i .
In conclusion, A flowchart describing the data analytical process employed in this study is available in Figure 1.

Clinical characteristics and survival analysis
We provide Table 1 to clinically depict the population that composes the cohort (n=112). We identify that most patients in the cohort are male (60%) with an average age of 64.7 years old. The average ICU stay was 38.3 days. Finally, 63.7% of all patients survived while the lowest reported mortality was in sepsis (30%) patients followed by sepsis COVID (32.4%) patients and septic shock patients (38%). We also report that the patients with septic shock have a higher count of neutrophils. Differences in the neutrophil counts were found to be associated with the severity of the disease; i.e., in the severe form of sepsis, the leukocytes count increased 1.62-fold and the increase in sepsis patients was 1.41-fold, while the leukocytes count remained more stable in septic shock (1.005-fold) patients. Next, we report the differential WBC count in sepsis, sepsis COVID, and septic shock (1.44, 1.15, and 1.08, respectively). Patients were assessed according to positive microbiological culture. We report that the COVID patients showed a smaller proportion of patients with viral and bacterial co-infection, i.e., superinfection (9 patients, 7 in LDMOF and 2 in HDMOF). The positive culture results for patients without COVID was found stable across both groups in sepsis and septic shock. Finally, we found the most common comorbidities to be hypertension, affecting 45% of the entire population of the cohort, followed by obesity (18%), and chronic obstructive pulmonary disease (15%). At the time of the study (i.e., September 2020 to March 2021), the circulating COVID-19 variant in the British Isles was B.1.1.7.
We also employed a survival analysis by days 28 and 90 ( Figure 2) to assess the mortality probability of sepsis, septic shock, and sepsis COVID patients. First, in Figure 2A, up to 28 days, the three groups of patients did not present significant differences in their survival probability (p=0.051); however, the low p value indicates a trend among the three groups' survivability rate, placing septic shock as the highest mortality rate after 28 days. When a 90-days analysis was considered ( Figure 2B), the survival rate among the patients showed statistical significance (p = 0.044), where the septic shock group presented the lowest survival probability and sepsis the highest.
Finally, to validate the binarization we performed, described in (18), we selected clinical parameters in our data that each, individually, represent the failure of a single organ. By following a logical expression that considers vasopressor as an exclusive (conjunction, AND) variable and platelets, bilirubin, creatinine, and P/F ratio as inclusive (disjunction, OR), our binarization resulted in an AUROC score of 0.91 followed by a Youden index of 0.71 (Supplementary Material S3).

Outcome of feature selection algorithms
To select a subset of biomarkers that explain the target variables (i.e., LDMOF and HDMOF), we performed a feature selection process. We chose three algorithms of distinct classes, namely Boruta, LR, and IG. Each application returned a different subset of biomarkers with varied lengths. We show in Figure 3, through panels A, B, and C, the outcomes of the Boruta, LR, and IG algorithms, respectively. The subsets of biomarkers obtained are as follows: i) PDL1, IL15, IL6, VCAM, IL1ra, IL1b, IL10, and CCL2 in Boruta; ii) MPO, VCAM, IL1b, VEGFC, IL17a, GMCSF, ANG2, CCL2, IL12, GRANB, and SGP130 in LR; and iii) PDL1, GRANB, IL15, ICAM, and IL1ra in IG.

Assessing the classification capacity of groups of biomarkers with different algorithms
We assessed the classification feasibility with five different machine learning algorithms, i.e., SVM, RF, CART, KNN, and Overview of the data analysis procedure employed in this study. In Figure 1, we show the three stages of the data analytical process employed in this study. First, 30 biomarkers from sepsis, sepsis COVID, and septic shock patients were obtained. Secondly, to reduce the number of variables, three algorithms are applied in the Feature Selection stage, i.e., Boruta, Lasso Regression, and Information Gain. The full outcome of the three algorithms is classified into an Artificial Neural Network (ANN). A second filter is applied to promote more reduction to the data using Exhaustive Search, whose outcomes are yet fed into ANNs in to compare their performance with the full outcomes of Boruta, Lasso Regression, and Information Gain. After running multiple ANNs, the prediction model is evaluated with SHAP.
ANN. We fed the classification models with the input variables identified by each one of our three FS algorithms ( Table 2). From that, we identified that the ANNs presented the most satisfactory predictability as its accuracy score outperformed the other methods.
The results displayed in Figure 4 indicate the detailed performance of the ANN models. To identify the best performing subset of biomarkers, we selected the model that presented the best error drop rate after 100 epochs, area under the curve (AUC), accuracy, precision, specificity, and recall. Next, the selected model was trained and tested with the subsets of biomarkers obtained in the FS stage. From that, the IG algorithm did not produce satisfactory results due to the imbalance of its prediction capacity, in which a low recall value was found (67.1%), indicating that it identified too many FNs in proportion of TPs, which is explained by the high error drop rate of the ANN model. Conversely, the ANNs trained/tested with the outcomes of LR and Boruta yielded satisfactory results, which is A B

FIGURE 2
Kaplan-Meier survival rate. Kaplan-Meier survival probabilities were identified. In (A), we show the Kaplan-Meier survival rate of sepsis, sepsis COVID, and septic shock patients after 28 days. Firstly, the data was identified as non-parametric (Shapiro-Wilk p = 1.377e-08) and the Kruskal-Wallis test was chosen to compare the averages (p = 0.051). In (B), we show the Kaplan-Meier curve for sepsis, sepsis COVID, and septic shock after 90 days, the data also follows a non-parametric distribution (Shapiro-Wilk = 6.036e-11) and the same Kruskal-Wallis test was employed to compare the averages (p = 0.044).

Model interpretation
To interpret the decision pattern of each ANN model, a SHAP approach was applied to the training dataset. Since the ANN trained/ tested with the biomarkers obtained by the Information Gain approach did not produce satisfactory classification metrics, we opted not to interpret this erroneous classification model. For both Boruta and Lasso regression, we analyzed the SHAP value for each input biomarker. We targeted the LDMOF class out of our train dataset and checked the contribution of each biomarker in predicting the class.
First, in Figures 5A, B, we show that both PDL1 and MPO were the biomarkers that mostly contributed to predictions in their models.
Next, in Figures 5C, D, we demonstrate the concentration of each biomarker in assigning the LDMOF label for each patient. From that, we see that there is a clear division between patients with high/low concentrations of input biomarkers (blue and red dots, representing each patient of the train set) getting negative and positive SHAP values, which directly affects the label assignment. To provide individual explanations for each biomarker, we targeted the ones that are positively correlated with LDMOF (i.e., VCAM, IL1ra, and IL4) and the biomarkers that are negatively correlated with HDMOF (i.e., PDL1, MPO, IL17a, and VEGFC). The remaining biomarkers did not have a clear separation between our target variables. Additionally, our two ANN models did not produce the same results regarding IL1b since it had different behaviors in each model. Therefore, the XAI approach of our tool enabled us to locate biomarkers with pro-and antiinflammatory activities. Biomarker selection using three feature selection algorithms. (A) indicates the feature selection process using the Boruta algorithm. The first set of obtained variables was submitted to a tentative fix method to deliver a more reliable subset. In the plot, the columns shown in green are the ones confirmed by the algorithm to be statistically significant and have higher importance in describing the data's label. (B) shows the results obtained by the feature selection using Lasso Regression. The x-axis of the figure indicates the log of lambda. Since Lasso Regression might be used as a classification model, the y-axis shows the AUC when including the number of variables shown at the top of the plot. The table below shows the variables stored in the lasso$lambda.min object, which corresponds to the variables with the best lambda values. The more distant a value is from zero, the more relevant it is for the predictor once that Lasso Regression sets the lambda = 0 to unimportant variables. (C) conveys the information gain derived from dataset entropy reduction achieved by the Information Gain algorithm. The vertical dashed line represents a threshold that considered the five most impactful biomarkers.

Conserved programmed death Ligand-1 and myeloperoxidase signals across patients with distinct manifestations of organ failure syndromes
We selected the two biomarkers flagged by the ANN (i.e., MPO and PDL1) to look for statistical differences between biomarker concentration and i) the three diseases of our cohort (sepsis, septic shock, and sepsis COVID-19) and ii) the two severity levels our binarization considered (LDMOF and HDMOF) ( Figure 6). Statistical significance (p = 0.007) was found only in the distinction of the MPO concentration among HDMOF and LDMOF patients ( Figure 6A), while the severity groups could not be statistically explained in terms of their PDL1 concentration ( Figure 6B). We only found statistical difference (p = 0.045) in using MPO as a distinguisher of sepsis and sepsis COVID ( Figure 6C), while the PDL1 concentration could not statistically differentiate sepsis, sepsis COVID, and septic shock ( Figure 6D).

Discussion
In our study, we could explain the severity of sepsis, sepsis COVID, and septic shock patients as a function of an unbalanced concentration of biomarkers. The input parameters of our function are subsets of biomarkers that explain a dysregulated host immune response. Moreover, we could isolate both MPO and PDL1 as the key contributors to the function.
The output parameters of our function are high and low degree of multi-organ failure (i.e., severity) of subgroups of patients that all meet sepsis-3 criteria, clinically placing the patients into a wider group, converging with past evidence (4). Our findings allowed us to immunologically place the patients from different disease models together as we failed to find major statistical significances in the concentration of PDL1 and MPO (Figures 6C, D) that distinguish sepsis, sepsis COVID, and septic shock all together. The only differences we could observe were in the MPO concentration between sepsis COVID and sepsis patients, and we argue that this is a result of the generally higher neutrophil count in bacterial induced sepsis (22). Finally, the SOFA score-based binarization we achieved is on par with results previously reported (18).
Next we showed that the WBC count increased with severity. Neutrophils, the most abundant WBC, were higher in the severe manifestations of the diseases we assessed. Moreover, the WBC and neutrophil count was lower in COVID sepsis, matching previous references of neutrophils being one of the most responsive cells toward bacterial infection (22). We argue that the lower occurrence of bacterial infection found in COVID (i.e., superinfection) patients contributed to their lower count of WBC and neutrophils.
Nonetheless, neutrophils remained an important immune cell to express cytokines, which might explain their response toward infection. Furthermore, no records of comorbidities influencing severity were found, except from cancer; we found the proportion of patients with cancer was higher in the severe form of the three disease models we assayed. The two key biomarkers we found (i.e., MPO and PDL1) are parts of important pathways in cancer immunotherapies (23) providing an opportunity for future studies of data science approaches for biomarkers involved in cancer.
Statistically, no significant differences were found that distinguished our patients based on their severity ( Figures 6A, B). It was previously reported that when statistics do not reach a satisfactory classification performance, Machine Learning (ML) might be a valid approach (24). We tried different ML approaches to classify our data. From five algorithms tested, only ANNs yielded a satisfactory classification. In fact, the robustness of this method was previously reported (24) as a solid way to find patterns in tabular data, among others. Additionally, the appropriate selection of the input variables is a key process in obtaining satisfactory results (25). In many cases, classification and regression models derived from lower-dimensional datasets benefit the downstream decision-making process (26,27). We further address this discussion by linking the appropriate selection of biomarkers with easily interpretable results in an information curation step. In this study, we were able to reduce a total of 30 biomarkers using three FS algorithms into subsets that conveyed satisfactory classification results in two instances. The three algorithms we selected belong to diverse classes of FS methods. In fact, they all have successfully been applied in reducing the complexity of ML inputs (20,(28)(29)(30)(31)(32) A B D C FIGURE 6 Statistical tests employed for PDL1 and MPO as characterizers of severity and disease model. Statistical significance tests comparing the averages of different groups of patients. In (A, B), we show the mean comparison of MPO and PDL1 (respectively) in distinguishing each one of the target labels (i.e., HDMOF and LDMOF). In (C, D), we employ the mean concentration MPO and PDL1 (respectively) in distinguishing sepsis, sepsis COVID, and septic shock.
without creating synthetic datasets. The results we obtained highlight the reduced complexity in employing datapreprocessing (DP) techniques. In fact, DP accounts for most of the workload involved in ML applications (33). We link the lack of success of the information gain algorithm in producing satisfactory predictability to the fact that this algorithm will only look for the association between input variables with a label (i.e., biomarkers and LDMOF/HDMOF). The other two FS algorithms will fit a model to determine the importance of each individual variable in predicting a label, granting them mathematical robustness. Therefore, we argue that a systematic evaluation of DP techniques such as the one here proposed is highly beneficial for developing insilico models.
After selecting optimal input biomarkers, we applied them in a classification system that uses this immunological information to predict the severity trajectory of critically ill patients. For that, our model, when fed with distinct subsets of biomarkers, could predict patients' severity. In fact, biomarkers have themselves been proposed as good predictors of a plethora of medical conditions (34)(35)(36); on the immunological side, they have been associated with pro-and anti-inflammatory responses (37) and we gathered evidence in support of our model succeeding to capture this. There are hundreds of biomarkers containing valuable information about the organic systems of the body and their functions. For their capacity to be interpreted, and consequently aid decision making and drug development, we argue that biomarkers related to a specific condition be systematically selected as we have proposed in this study.
We interpret our ANN model by linking the concentration of biomarkers in explaining LDMOF. First, Multiple Organ Dysfunction and even death have been reported to be associated with increased levels of VCAM-1 in adults and neonates diagnosed with sepsis (38)(39)(40). Similarly, our model shows a negative correlation between increased VCAM levels and LDMOF development. Next, IL1ra was used in a recombinant treatment and was successful in reducing levels of mortality (41). Our last biomarker positively associated with LDMOF is IL4, which was reported to act together with IL6 to induce Th2 cells and macrophage differentiation (42). Lower concentrations of this biomarker were associated with lower mortality of severe sepsis patients (43).
We spotted four biomarkers that are negatively correlated with LDMOF, granting them a negative effect on a patient's favorable outcome. First, low concentrations of IL17a have been reported as a good predictor for mortality in sepsis caused by distinct pathogens (44,45). Next, septic shock and sepsis can lead to hypoxia due to tissue hypoperfusion (1). A transcription factor named hypoxiainducible factor 1-alpha (HIF-1a) accumulates in cells under hypoxic conditions and can upregulate the expression of VEGF and PD-L1 (46,47). The activated HIF pathway can also trigger the activation of innate immune cells, including macrophages, dendritic cells, neutrophils, and natural killer cells (48). Furthermore, increased levels of PD-L1 will suppress the adaptive immune response by inhibiting the proliferation and activity of the CD4+ effector T cells and enhancing the differentiation of Tregs (49, 50). This persistent inflammatory condition caused by innate immune activation along with suppressed adaptive immune response contributes to hypoxia-induced organ failure. It has also been reported that PD-L1 knock-out animals show a better survival rate after a septic challenge compared to wild-type animals (51). VEGF is also known to induce VCAM-1 expression, and its upregulation is correlated with damaged vascular endothelium and organ dysfunction (40,52). Thus high levels of PD-L1 and VEGF can be considered key markers of multi-organ failure. Finally, myeloperoxidase (MPO), an enzyme produced by neutrophils, is found to be increased in severe patients suffering from septic shock, despite having no significant difference in neutrophil count (53). In another study, MPO levels during the early stages of sepsis were found to be negatively correlated with patient survival (54).
The mapping of the XAI model provided meaningful information on the course of dysregulated immune responses and also converged with clinical interpretations regarding the neutrophil count of the disease models (i.e., neutrophils tend to increase with severity). A compelling example is the proinflammatory function our model attributed to MPO. This enzyme is primarily produced by the granulocytes of neutrophils. The overexpression of MPO generates harmful chemicals that have a detrimental effect on organ inflammation (11). Our model also identified PDL1 as a pro-inflammatory protein. The blockade of the binding between PDL1 and PD1 might inhibit lymphocytes from apoptosis. The upregulation of PDL1 by neutrophils is increased in sepsis as the higher migration of these cells might allow them to be trapped in the lung vasculature (55, 56). Therefore, with more neutrophils expressing PDL1, immunosuppressant effects start to occur with the death of neutrophils. The highest counts of neutrophils were found in the severe manifestations of the three disease models we investigated. Therefore, we can track the course of the multi-organ failure syndrome of our patients with increased neutrophils leading the overexpression of MPO and PDL1.
In-silico models are an efficient paradigm of experimentation. Compelling examples are found in big data-based applications that have been assisting in several areas in the medical sciences, such as predicting heart attacks (57), telediagnosis (58), and preventing disease outbreaks (59), among others. A guideline proposed by the Center for Drug Evaluation and Research (CDER) (60) shows that the development of therapeutics might initiate with screening characteristics that indicate biological processes. Here we propound to employ data science as an initial step for screening biomarkers, enabling the gathering of solid mathematical evidence for linking concentrations of biomarkers with patients' severity.

Conclusions
In the present study, we built a function whose input is the concentration of biomarkers and the output is the level of severity of a patient. For our goal to be achieved, a systematic data mining procedure enabled us to identify the upregulation of PDL1 and MPO as good predictors of severity in sepsis (viral and bacterial induced) and septic shock patients. After interpreting the results both clinically and immunologically, we found that there is solid medical and biological evidence for why the upregulation of PDL1 and MPO is a major driver of severity. To this extent, we posit that data mining routines such as the one we proposed be used to identify the biomarkers that can function as part of an early diagnosis system.

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 in the article/Supplementary Material.

Ethics statement
The studies involving human participants were reviewed and approved by Health Research Consent Declaration Committee (HRCDC) under the register REC: 2020-05 List 17 and project ID 0428. The patients/participants provided their written informed consent to participate in this study.  . The funders of the study had no role in the study design, data collection, analysis, or interpretation, writing of the report, or decision to submit for publication. DK is the Canada Research Chair in Translational Vaccinology and Inflammation.