No clinical trials have been conducted for this study. The data used in this study were obtained from the CoV-2-STs-001 study, assessing the safety and efficacy of CoV-2-STs compared with the SoC in hospitalized adults with severe COVID-19 during delta variant predominance [12]. In this trial, hospitalized patients with high-risk of COVID-19, within the first 6 days of symptom onset, with pneumonia or oxygen saturation on air, equal to or less than 94%, CD3+ cells less than 650 cells/μL, and at least 1 elevated serum biomarker (C-reactive protein [CRP], lactate dehydrogenase [LDH], ferritin, and D-Dimers) were randomly assigned (2:1) to receive either CoV-2-STs and SoC (CoV-2-STs-group; arm A: n=57 patients) or SoC-alone (SoC-group; arm B: n=30) and were followed for 60 days. Both treatment arms received the same SoC. The same baseline SoC, which comprised dexamethasone and remdesivir as part of the recommended therapy for hospitalized patients with severe COVID-19, was given to both the CoV-2-ST+SoC and SoC-only groups. Comparability was maintained as the investigation found no appreciable variations in how these therapies were administered between the 2 groups. However, depending on clinical assessment, there were some differences in the usage of extra supportive therapies. At the clinician’s discretion, several patients in the SoC group received tocilizumab when their cytokine levels continued to be high.

The cells were infused within a maximum of 24 h after randomization (day 0 for both arms). Clinical data were assessed and recorded daily. Laboratory assessments, including (1) levels of CRP, LDH, ferritin, D-Dimers, and white blood cells (WBCs); (2) number of circulating CD3+, CD4+, CD8+, and CD56+ cells; and (3) CoV-2-STs, SARS-CoV-2–specific neutralizing antibodies (SARS-CοV-2-NAbs), and serum cytokine measurements, were performed on days 0, 3, 5 (±2), 10 (±2), 15 (+3), 22 (+3), 30 (+3), 45 (+3), and 60 (+7).

Serum cytokines (interleukin [IL] 6) and tumor necrosis factor-α (TNF-α) were tested to monitor systemic inflammation in patients. The frequency of CoV-2-STs was measured using interferon-gamma (IFN-γ) enzyme-linked immunospot (ELISpot) assays, which detect IFN-γ production after antigen stimulation. Direct ex vivo ELISpot testing showed T cells’ functional responsiveness to SARS-CoV-2 antigens.

For the scope of this study, only the biomarker values assessed on days 0 and 5, and their correlation with the outcomes of recovery (defined as a patient score ≤3 in the 8-category World Health Organization [WHO] ordinal scale) and survival were considered [25]. The patient population was subjected to an “as treated” analysis, including in arm B, 1 patient who was initially randomized to arm A but was intubated before receiving the CοV-2-STs (arm A, CoV-2-STs+SoC n=57 patients, and arm B, SoC-group n=30 patients).

The primary objective of this research endeavor was to develop, by using robust data analysis techniques, a medical tool capable of predicting with high accuracy the ultimate treatment outcomes of severe COVID-19. The treatment modalities assessed were the SοC alone (including dexamethasone and remdesivir) or the adoptive transfer of CοV-2-STS in addition to SοC, from the randomization of patients (day 0) until study completion (day 60).

Our strategy involved a multistep process outlined in Figure 1. Initially, we pinpointed a critical, early postenrollment time point (day 5) where significant biomarker differences were observed between the 2 study groups (arms A and B), thus implying that the treatment impact was visible by day 5 on patients’ biomarker profiles. These differences were detected by analyzing the biomarkers’ descriptive statistics at baseline and on day 5 and calculating the rate of change for each biomarker over the initial 5 days. We evaluated a panel of supervised learners—L2-regularized logistic regression, penalized or shrinkage linear discriminant analysis (LDA), random forest, and gradient boosting—to predict day-60 outcome from day-5 measurements. The end point was defined a priori as a binary composite: favorable (alive in remission) versus critical (deceased or alive-with-disease). Modeling was performed separately by treatment arm. All preprocessing occurred within cross-validation (CV) folds to prevent leakage. Numeric features were median-imputed and standardized (z scores); categorical features were mode-imputed and one-hot encoded. We used the same stratified CV folds for all models with nested hyperparameter tuning [26,27].

To substantiate model choice, we compared 5 different models, including LDA solved via singular value decomposition (svd), shrinkage or regularized LDA, regularized logistic regression, random forest, and gradient boosting using identical stratified k-fold CV splits for all methods [28,29]. Within each training fold, we performed inner-fold tuning, and all preprocessing (imputation and standardization) was confined to the training data to prevent leakage. We reported out-of-fold (OOF) performance for each model, including area under the receiver operating characteristic curve (AUROC), area under the precision-recall curve (AUPRC), Brier score, sensitivity, specificity, and F₁-score.

Feature attribution was computed within CV folds and then aggregated across folds. For logistic regression, the standardized coefficients could be used. For shrinkage-LDA, the standardized discriminant loadings from linear discriminant 1 (and linear discriminant 2, where defined) could be reported as the primary importance measure. For tree-based models (random forest and gradient boosting), model-agnostic permutation importance and Shapley Additive Explanations values using fold-fitted models could be computed [30]. In our case, based on the best model selected, the cross-model concordance of top features is summarized, and univariate confirmation (Kruskal-Wallis [KW] for continuous variables; chi-square for categorical) is provided to contextualize multivariable findings [31,32]. Using the above workflow, we trained 2 arm-specific models: ML(CoV-2-STs+SoC) for arm A and ML(SoC-only) for arm B, each using day-5 predictors. Comparative performance summaries and calibration diagnostics are presented in the Results section.

The ultimate goal of our efforts was the development of a computational tool capable of predicting a patient’s most probable outcome based on their initial characteristics and the treatment applied, thus enabling optimal therapeutic decision-making. Using this tool, we conducted simulations to assess the potential benefits of administering CoV-2-ST treatment to patients who initially received only the SoC (arm B) and the potential adverse outcome of treating with SoC only to patients who had initially received CoV-2-STs. For these simulations, we applied a Monte Carlo method to estimate day-5 biomarker values, considering the distribution of biomarker actual rates of change from day 0 to day 5, in each group. This approach allowed us to generate 1000 hypothetical scenarios for each group and predict clinical outcomes at day 60 using our ML(SoC-only) and ML(CoV-2-STs+SoC) models, respectively. Our approach for both ML models predicting patients’ outcomes with CoV-2-STs+SoC or SoC-only treatment is summarized in Figure 1.

The input data for each group encompassed demographic, clinical, and laboratory parameters, including age, gender, Karnofsky score (KS), comorbidities, serum levels of IFN-γ, IL-10, IL-2, IL-6, TNF-α, SARS-CοV-2-NAbs, circulating CD3+, CD4+, CD56+, CD8+, and CoV-2-STs, and levels of CPK, CRP, D-Dimers, ferritin, LDH, and WBC. The dataset was subjected to standardization using the Standard Scaler. The correlations between the features were quantified (Figure S1 in Multimedia Appendix 1). However, highly correlated features were not removed because LDA classifications transform the input variables into a set of orthogonal uncorrelated variables [33]. Also, it is known that collinearity does not impact the performance of the LDA classification. Both SoC and CoV-2-STs groups, k-fold CV with k=5 was used, ensuring robust evaluation by partitioning the dataset into 5 subsets and iteratively using 4 folds for training and 1 for validation in each iteration. For the classification task, the classifier was used based on patients’ day-5 biomarker levels, along with baseline characteristics, such as age, comorbidities, and KS, and used to categorize the patients into distinct outcomes, such as (1) alive and well, (2) deceased, and (3) alive with disease. Among disease severity variables (KS, oxygen supply, and 8-point category WHO Ordinal Scale), KS was selected to be incorporated in the classification model, as all variables were highly correlated [12,34]. The model performance was assessed using the following metrics: AUPRC, AUROC, sensitivity, specificity, F₁-score, and Brier score. Because the outcome is imbalanced, we prespecified AUPRC as the primary discrimination metric. We also report AUROC as a threshold-free complement. To summarize performance at an actionable operating point, we compute sensitivity, specificity, precision or recall, and F₁-score at the fold-wise threshold that maximizes F₁-score on the validation split. Calibration is assessed with the Brier score. We summarize metrics as mean (SD) across the same stratified CV folds. Accuracy is not used as a headline metric because it can be inflated by the majority class in imbalanced datasets.

For categorical variables, prevalence and percentage values were used to calculate descriptive statistics. Medians and IQR were used for continuous variables, and median and range were used for variables with skewed distributions. For nonparametric comparisons, the KW and Mann-Whitney U tests were used [31]. The Fisher exact test or the chi-square test, as appropriate, was used to test categorical comparisons [34].

The data used in this study were obtained from the CoV-2-STs-001 study (Clinical Trials.gov NCT05447013, EudraCT:2021-001022-22), which was conducted in accordance with the Declaration of Helsinki. Approval for the clinical trial was obtained from the National Organization for Medicines and jurisdictional ethics committees (the National Ethics Committee and the institutional review boards of George Papanikolaou Hospital and Hippokrateion Hospital; IRB 5th/03-06-2021) after thorough discussions and protocol adjustments on eligibility criteria and safety assessments and an interim and final review of the phase 1 study data (National Ethics Committee) before beginning the phase 2 study. Continuous safety surveillance by pharmacovigilance was also integrated into all study phases.

Written informed consent was obtained from all participants before enrollment, including consent for the use of deidentified clinical and biomarker data for subsequent research analyses. All data used in this post hoc ML analysis were fully anonymized before analysis, and no personally identifiable information was accessible to the investigators performing the computational analyses, thereby ensuring participant privacy and confidentiality. Participants did not receive financial compensation for participation in the present secondary analysis.

The baseline demographic and clinical characteristics of study participants were well-balanced between the 2 groups and are summarized in Table 1. The corresponding intention-to-treat analysis has been reported in detail elsewhere [12]. Importantly, while at baseline there were no differences in the levels of tested biomarkers between patients treated with CoV-2-STs+SoC and SoC, apart from CPK (P=.01; Figure 2A), on day 5 post infusion, the rate of changes in the levels of CD3+ T cells, natural killer cells (CD56+ natural killer [NK] cells), CD8+ T cells, and CoV-2-STs and from baseline, were significantly different between the 2 arms (Figure 2B), with increases favoring the CoV-2-ST group (Mann-Whitney P values of .04, <.001, .01, and .02, respectively).

The distinctive immunological responses elicited by CoV-2-STs as opposed to SoC treatment on day 5 identified this day as a critical time point. On day 5, CoV-2-STs selectively modulated certain biological markers, leading to enhanced adaptive immunity (as evidenced by the increases in T cells and CoV-2-STs) while preserving innate immunity (indicated by maintained levels of NK cells). The day-5 data subsequently played an instrumental role in the development of our ML model.

To gain deeper insights into how the 2 treatments may differentially modulate the immune response and the biomarker components interact, we analyzed the correlations among all biomarkers on day 5 post treatment for each group (Figure S1 in Multimedia Appendix 1). As expected, a significant correlation (exceeding 0.7) between CD3+ and CD4+ cells was observed in both groups, suggesting a predominant presence of helper T cells in circulation. This aligns with the established critical role of CD4+ T cells in immune regulation and viral clearance in COVID-19. In group A, the observed correlation between IFN-γ and IL-2 indicates a robust T cell response. IL-2 is essential for T cell proliferation, while IFN-γ is a crucial cytokine produced by activated T cells. This correlation implies that the administration of CoV-2-STs enhances T cell–mediated immune defense against the virus. Conversely, in group B, significant correlations of proinflammatory cytokines IFN-γ and IL-6 suggest an innate and broad inflammatory reaction aiming at combating the virus and resolving acute inflammation. Additionally, the correlations between proinflammatory and anti-inflammatory cytokines, such as IFN-γ and IL-10 or IL-6 and IL-10, may represent mechanisms to mitigate potential collateral damage from an intense immune response. Furthermore, a strong correlation between CPK and TNF-α was observed, linking muscle impairment in severe COVID-19, such as myalgia, rhabdomyolysis, or myositis, to a stress response involving proinflammatory TNF-α [35-37]. Overall, these findings underscore qualitatively distinct immune modulations by the CoV-2-ST treatment compared with SoC.

Our objective was to use early posttreatment biomarkers from day 5 to predict day 60 clinical outcome and highlight markers most informative for prognosis. Preprocessing included fold-wise transformation and standardization to avoid leakage. Because of severe class imbalance and to stabilize estimation, we defined a binary composite end point with favorable (alive and healthy) versus critical (deceased or alive-with-disease). Models were fit with stratified CV and evaluated separately within each cohort or treatment arm.

We first benchmarked regularized logistic regression, shrinkage-LDA, random forest, and gradient boosting using the same stratified CV folds with nested tuning. We report AUROC, AUPRC (primary), sensitivity, specificity, F₁-score, and Brier, with CIs from repeated CV and calibration curves. The index model was chosen as the classifier with the highest mean AUPRC and competitive calibration (Brier). Feature importance was then quantified on the index model using the method appropriate to its class—standardized coefficients for regularized logistic regression; signed, standardized discriminant loadings for shrinkage-LDA; and Shapley Additive Explanations with in-fold permutation importance for tree-based models—with cross-model concordance. On arm A (Table 2), shrinkage-LDA performs competitively with the best nonlinear alternative. Its AUROC (mean 0.874, SD 0.052) is within fold-to-fold variability of random forest (mean 0.881, SD 0.079), with AUPRC and F₁-score similarly close (0.832 and 0.842, respectively), and a respectable Brier score (0.137). Thus, shrinkage-LDA is not dominated by the more complex models while retaining a transparent linear decision rule. On arm B (Table 3), random forest attains the highest mean discrimination and F₁-score (mean AUROC 0.739, SD 0.042; mean AUPRC 0.862, SD 0.077; mean F₁-score 0.836, SD 0.048), whereas shrinkage-LDA delivers the highest sensitivity (mean 0.950, SD 0.112) at the cost of lower specificity and a higher Brier (mean 0.294, SD 0.101). Given the small-n setting and class imbalance, this sensitivity-favoring profile may be desirable when false negatives are clinically costly, and we therefore report shrinkage-LDA alongside the index model. Because false negatives are costly, the best model was specified as the one that achieves noninferior discrimination and calibration (within CV variability) and the highest sensitivity at a clinically relevant threshold. Shrinkage-LDA satisfies this on arm A and yields the highest sensitivity on arm B, coupled with its interpretability and regularization. Shrinkage-LDA was selected as the primary clinical model and reports its signed, standardized loadings as the main feature-importance summary. KS was chosen as a representative variable over other disease severity variables, as mentioned earlier. While random forest offered superior overall discrimination (AUROC 0.739), we prioritized shrinkage-LDA for the clinical decision tool because, in the context of severe COVID-19, the “cost” of a false negative (missing a patient who will die) is far higher than a false positive (overtreating a patient who would have recovered). LDA’s near-perfect sensitivity (0.95) ensures patients with high risk are captured, even at the expense of lower specificity. Factors were initially ranked in terms of importance for each group, followed by further filtering to consider only those that significantly differed between the 2 outcome groups (alive and well, and critical) by KW test (Table S1 in Multimedia Appendix 1). According to this analysis, key determinants for the day-60 outcome in the CoV-2-STs+SoC arm included (1) the KS; (2) the levels of CD3+, CD4+, CD8+, and CD56+ cells; (3) the levels of circulating CoV-2-STs; and (4) the levels of IL-6, LDH, ferritin, CRP, and WBC, while in the SoC-arm only the circulating CoV-2-STs, CD3+, CD56, and the CRP levels were identified as significant (Table S1, Figures S2 and S3 in Multimedia Appendix 1). Based on these data, we proceeded to further analyses to define suggested cutoffs for critical variables, presented in Table 4, according to the corresponding treatment and outcome.

Given the limited sample size and class imbalance, we modeled a binary end point—favorable (alive and well) versus critical (deceased or alive with disease)—using repeated stratified 5-fold CV (10 repeats) with identical splits across strategies and all preprocessing confined to the training fold. We evaluated shrinkage-LDA under 3 data-retaining imbalance strategies: empirical priors (no resampling), balanced priors (priors=0.5-0.5; LDA analog of class weighting), and Synthetic Minority Oversampling Technique (SMOTE) after one-hot (pragmatic oversampling). In each fold, the F₁-optimal threshold was chosen on the validation split; training or validation confusion matrices were computed and aggregated. Because of imbalance, AUPRC was prespecified as the primary metric; we also report AUROC, sensitivity, specificity, precision, F₁-score, and Brier, with 95% bootstrap CIs (B=2000) computed on OOF pairs while refitting the threshold within each bootstrap. Arm-specific OOF summaries, confusion matrices, and fold means are in Tables S2-S4 in Multimedia Appendix 1 (arm A) and Tables S5-S7 in Multimedia Appendix 1 (arm B).

In the CoV-2-STs+SoC arm (n=57), OOF discrimination was similar across modeling strategies (Table S2 in Multimedia Appendix 1): Brier=0.14‐0.16, AUROC=0.860‐0.880 (overlapping 95% CIs), and AUPRC=0.740‐0.780. At the F₁-optimized operating point (Table S3 in Multimedia Appendix 1), combined validation rates were strong and well-balanced (F₁-score=0.810‐0.840 with precision=0.740‐0.790, specificity=0.830‐0.870, and sensitivity=0.890‐0.910). Fold-level summaries (means and SD; Table S4 in Multimedia Appendix 1) recapitulated these patterns. Among priors, balanced priors tended to yield the highest sensitivity (approximately 0.910), whereas empirical priors offered slightly higher precision and F₁-scores (approximately 0.83‐0.84). SMOTE after One-Hot Encoding performed comparably to these settings and did not provide a clear advantage.

In the SoC arm (n=30), OOF performance again showed overlapping CIs across strategies (Table S5 in Multimedia Appendix 1), with AUPRC=0.78‐0.81, AUROC=0.72‐0.76, and Brier=0.22‐0.27. Aggregated validation confusion (Table S6 in Multimedia Appendix 1) emphasized very high sensitivity (approximately 0.95 across strategies) at the cost of lower specificity (0.58‐0.62), resulting in precision=0.77‐0.79 and F₁-score=0.85‐0.87. Fold-level means (Table S7 in Multimedia Appendix 1) were consistent (AUPRC=0.86‐0.88, AUROC=0.72‐0.75, Brier=0.25‐0.27, and F₁-score=0.86‐0.87). Balanced and empirical priors were nearly indistinguishable, and SMOTE after One-Hot Encoding again produced no material gain—an overall profile that is expected given the smaller, more imbalanced sample.

These quantitative results align with the OOF visualizations in Figure 3. Calibration plots (A and B) show predicted probabilities tracking observed risk, while receiver operating characteristic (C and D) and precision-recall curves (E and F) indicate good separation in arm A and fair-to-good separation in arm B. Decision-curve analyses (G and H) demonstrate positive net benefit for model-guided decisions over “treat-none” and over “treat-all” across clinically relevant threshold ranges, supporting the practical value of risk-based treatment selection in both arms.

After fitting the 2 arm-specific models, LDA-ML_{CoV-2-STs+SoC} and LDA-ML_SoC, we derived a predictive “what-if” tool to estimate how each patient’s predicted day-60 outcome might change under the alternative treatment (Figure 4). This analysis is associational and does not identify causal effects. Also, this tool will enable the identification of patients who have low probabilities to recover with SoC-only, while representing the best candidates for adoptive immunotherapy with CoV-2-STs or other alternative treatments. As an input, we used the baseline characteristics and biomarker values of 21 predictors in total (age, gender, KS, comorbidities, IFN-γ, IL-10, IL-2, IL-6, TNF-α, CD3, CD4, CD8, CD56, SARS-CoV-2-NAbs, CoV-2-STs, CPK, CRP, D-Dimers, ferritin, LDH, and WBC). Day 5 was prespecified as the early posttreatment anchor—the earliest uniform assessment across arms following therapy administration—because it captures the initial on-treatment biomarker dynamics and is the time point used for all prognostic models. Using baseline (day 0) and this posttreatment day-5 time point, we first computed the vector of 5-day changes (Δ=day 5–day 0) across biomarkers.

For the arm-swap Monte Carlo scenario analysis, day-5 biomarker values were not taken directly from the observed data but were generated from 4 variance models for the 5-day change rates: (1) additive independent and identically distributed (IID) bootstrap, (2) additive IID normal (clipped to the donor 1%‐99% range), (3) additive multivariate normal with a Ledoit-Wolf shrinkage covariance, and (4) multiplicative IID bootstrap on percent-per-day changes. Specifically, for each arm, we simulated day 5 values as (day 5_sim=day 0+Δ̃) where Δ̃ was drawn from the corresponding variance model for that arm. Simulated day-5 values were clipped to physiologic ranges when necessary. For each individual, we generated n=1000 Monte Carlo scenarios per arm. Each simulated day-5 profile was then scored with the donor arm’s LDA model, yielding a predicted probability of the critical class at the day-60 outcome, and for each patient, we stored the mean predicted probability across Monte Carlo simulations.

Swap summaries are evaluated at a single fixed threshold of 0.50. A 0.5 probability threshold is the standard “neutral” cutoff for binary classifiers and corresponds to a decision threshold where false positives and false negatives are weighted equally in simple cost terms. Using a constant 0.5 avoids threshold-tuning bias, sidesteps calibration differences between donor and recipient arms, and keeps the reported critical class easy to interpret. We provide baseline-risk quartile subgroup counts, and repeat the swap with CV-trained donor models to obtain mean and SD class-1 counts across the Monte Carlo simulations. Swapped-arm predictions are predictive scenario analyses, not causal effects.

When arm A receives SoC-like early dynamics (Figure 5A), at a stringent 0.5 cutoff, the predicted number in the critical class depends on the variance model but is always ≥ the observed 35.1%. With the multivariate normal (MVN) model, the model calls 54 (SD 3; across CV refits) out of 57 patients, corresponding to 94.7% of the cohort and approximately 34 additional patients relative to the observed outcome. With IID additive generators the increase is moderate, 31/57 (SD 6; 54.4% of the cohort; approximately 11 additional critical patients), and additive IID normal, 31 (SD 5) critical cases out of 57 patients (54.4% of the cohort; approximately 11 patients extra). With multiplicative IID, the estimate is 21 (SD 6) out of 57 (36.8% of the cohort; approximately 1 patient extra). SDs across 5-fold CV refits are small relative to the between-model gaps, indicating that modeling assumptions, not refitting noise, drive the shift at this threshold.

When arm B receives CoV-2-STs+SoC early dynamics, at a stringent 0.5 cutoff, the results are highly sensitive to correlation (Figure 5B). Under MVN, the model calls 30 (SD 1; across CV refits) out of 30 patients, corresponding to 100% of the cohort, and approximately 12 patients extra. Under IID additive assumptions, the predictions drop well below trial, 11 (SD 2) out of 30 patients (36.7% of the cohort; 7 patients less), and additive IID normal, 11 (SD 3) out of 30 patients (36.7% of the cohort; 7 patients less). With multiplicative IID, the estimate is 9 (SD 2) out of 30 patients (30% of the cohort; 9 patients less). Again, the SDs are small compared with the MVN-IID separation, underscoring that retaining cross-marker comovement is the key determinant at high thresholds. The divergence of the MVN model in arm B (predicting 100% critical outcomes) is likely attributable to the instability of estimating a high-dimensional covariance matrix from the small sample size (n=30), which can amplify noise. Consequently, the independent (IID) bootstrap models likely provide a more reliable estimate for this specific small-cohort scenario.

These are predictive scenario analyses, not causal effects. They assume harmonized day-0 measurements, transportability of the donor from day 0 to day 5, change distributions to recipients, and no interference between patients; thresholds are learned on the donor arm, so cross-arm counts should not be read as treatment effects. The sensitivity panels across the 4 variance models are therefore the right way to present robustness (from CoV-2-STs+Soc to SoC-only) and dependency (from SoC-only to CoV-2-STs+SoC) without making causal claims.

We need to acknowledge here that cellular immune responses are not routinely assessed in health care facilities. Therefore, we developed an alternative version of the tool to facilitate its broader applicability, which excludes the circulating CoV-2-STs. This model still maintains the same characteristics as the model including all biomarkers (Figure S4 and S5 in Multimedia Appendix 1).

In previous studies, both our group and others have demonstrated a strong association between CoV-2-ST responses and favorable outcomes [12,38,39]. Furthermore, in a recent study, we showed that incorporating CoV-2-STs generated from recovered donors into the standard care regimen—including remdesivir and dexamethasone—markedly improved recovery and survival rates in patients with severe COVID-19, caused by the delta variant, compared with the SoC alone. This method of adoptive immunotherapy led to a rapid increase of circulating CoV-2-STs in vivo and contributed to the overall rebuilding of the immune system. Expanding on these findings, we here used ML for a post hoc analysis of the already-published clinical trial data to identify pivotal time points, characteristics, and biomarkers that are decisive for the outcome in severe COVID-19 cases under both treatment approaches [12]. Our analysis pinpointed day 5 post enrollment as a critical juncture at which distinct immune responses to CoV-2-ST administration compared with SoC become apparent, especially in the kinetics of T lymphocytes, circulating CoV-2-STs, cytotoxic T lymphocytes, and NK cells. This early distinction indicates a superior immune modulation by CoV-2-STs, impacting both innate and adaptive immunity, and corresponds with our clinical findings of faster recovery from lymphopenia, maintenance of innate immunity homeostasis, and enhancement of SARS-CoV-2–specific immunity in patients treated with CoV-2-STs in addition to SoC [12].

Furthermore, by using regularized or shrinkage LDA, we determined key predictors of patient outcomes by day 60, underlining the role of SARS-CoV-2–specific immunity in managing severe COVID-19 cases. Our LDA models, tailored for forecasting responses to CoV-2-STs+SoC or SoC-alone treatment, facilitated the development of a computational tool aimed at guiding therapeutic decisions. This tool, by considering initial patient characteristics and biomarker levels, seeks to identify individuals among those with severe COVID-19 who may not respond to SoC favorably, thus becoming prime candidates for alternative interventions, such as adoptive immunotherapy with CoV-2-STs. By simulating over 1000 potential disease progression profiles and using our models to compute survival probabilities, we present a new method for choosing treatment strategies that optimize patient outcomes. Testing this tool in hypothetical scenarios, we have further validated its potential for practical application, showcasing its adaptability to treatment modifications and its predictive accuracy across diverse patient responses, thereby highlighting its value in personalized medicine for severe COVID-19.

Our study benefits from the inclusion of a consistent sample of delta variant, high-risk individuals, mirroring real-world scenarios. Additionally, the availability of longitudinal, clinical, and laboratory data prospectively collected throughout our previous phase I and II clinical trial adds strength to our analysis, despite its post hoc nature. Leveraging these data, our prognostic approach enabled us to make early and precise estimates of mortality upon admission. Although reported ML models have provided promising prognostic implications, they are still hampered by certain limitations, most commonly the lack of capability to predict the outcome at, or early after, admission. Our rapid risk assessment tool addresses this gap, empowering health care professionals to identify patients who are highly susceptible to experiencing unfavorable outcomes under routine treatment alone. This allows for immediate and individualized interventions, such as improved care protocols and monitoring practices, and the consideration of administering CoV-2-STs. Notably, our research pioneers the integration of circulating CoV-2-STs into an ML-based predictive tool, enhancing the precision and applicability of our prognostic model. This approach aligns with the learning health system principles successfully implemented in large-scale COVID-19 surveillance registries [40].

Our study also has certain limitations that need to be acknowledged. First, the patient sample size was limited, as the original clinical trial was not specifically tailored to accommodate an ML analysis, whereas the retrospective nature of our research precluded prospective evaluation in patient cohorts. Because this is a post hoc analysis of a trial not purpose-built for ML, our findings should be interpreted as predictive and exploratory rather than causal or deployment-ready. To strengthen generalizability, we outline external and prospective validation in independent cohorts with a locked model, prespecified thresholds, and ongoing calibration monitoring. Second, since the initial trial only included severe cases, our tool’s evaluation was confined to the severe COVID-19 context, limiting its broader applicability. Furthermore, vaccination status was not included as a variable in our computational tool, due to the low vaccination rates among study patients (10/57, 17.5% of the CoV-2-ST group and 8/30, 26.7% of the SoC group). Since significant protection was offered by vaccination against severe illness, the majority of the eligible patients with severe COVID-19 were unvaccinated. To prevent potential bias stemming from underrepresented vaccinated patients, we chose not to include vaccination status in the tool. Future work should explore the incorporation of interoperable data standards to enhance reproducibility, as exemplified in large European vaccine trials [41].

While directly useful, the ML models depend heavily on the quality and breadth of the data used. If the data are not representative of the broader population or if key variables are missing, the predictions may not be accurate or generalizable. Accordingly, performance may be sensitive to the data context and temporal drift, and transportability beyond severe COVID-19 should be verified with temporal or external validation. Additionally, treatment outcomes can be influenced by multifactorial and dynamic patient conditions not fully captured in the model. The ML COVID-19 studies mentioned in the Introduction section use a variety of data types (eg, imaging and clinical parameters) and are often designed to process and analyze complex datasets quickly, providing critical insights in fast-moving clinical scenarios like a pandemic [19,20,23,24,42]. Their use of advanced algorithms can detect patterns not readily apparent to humans, which is a significant advantage in diagnostic and prognostic contexts. While our treatment comparison model provides specific actionable data for clinical decisions between 2 treatments, the broader ML studies offer insights into disease characterization, risk stratification, and initial management strategies. Future work will include prospective, preregistered evaluation and harmonized data standards to improve reproducibility and clinical uptake. A comprehensive approach in clinical settings might involve integrating insights from both types of models—using ML-driven diagnostics and risk assessments to initially guide treatment choices and then applying the treatment comparison model to select the optimal therapy based on clinical outcome predictions. In essence, while the referenced ML works provide broad diagnostic and prognostic tools useful in managing COVID-19, our model comparing treatments between CoV-2-STs and SoC offers targeted, actionable insights specifically designed to optimize therapeutic outcomes. Both have distinct roles in clinical decision-making, demonstrating the diverse applications of ML in health care.

This study focuses on the application of ML methods to COVID-19 clinical management, introducing several innovations. Previous models have successfully predicted mortality, severity, or intensive care needs based on early clinical parameters, imaging data, and laboratory results [20,23,24,42]. However, these models often rely on static baseline data and do not incorporate dynamic immune monitoring, particularly T cell–specific parameters. Most ML COVID-19 models focus on population-level predictions rather than individualized therapeutic decision-making, and few address which patients will benefit from emerging advanced therapies. The COVID-19 ML application focuses on clinical interpretability and direct applicability, using LDA instead of opaque prediction tools. LDA generates transparent decision boundaries and provides clear insights into feature contributions, making it suitable for structured biomedical datasets. Its reliance on linear combinations of input features and ability to identify discriminative variables make it suitable for structured biomedical datasets. This aligns with calls from the clinical ML community to prioritize explainable models that foster trust and adoption among health care practitioners [33,43,44]. The work combines ML, immunological profiling, and outcome simulation to enhance precision medicine in severe COVID-19.

Future research possibilities will prioritize ensemble strategies to improve the robustness and clinical utility of our COVID-19 prognostic models. First, we will evaluate bagging with gradient-boosted base learners alongside principled preprocessing (eg, outlier detection and information-gain feature selection) to handle small, noisy clinical cohorts [45]. In parallel, we will develop hybrid ensembles (stacking or soft-voting) that combine linear, tree-based, and margin classifiers to capture complementary signals across biomarkers [46]. All candidates will be benchmarked under an identical stratified CV with AUROC or AUPRC, Brier, calibration (pre-post isotonic), decision curves, and thresholded confusion matrices, followed by temporal or external validation to test transportability and drift. Finally, we will pair ensembles with clinician-oriented explainability (global or local summaries and reliability tables) and outline a prospective evaluation to assess workflow impact.

In conclusion, we provide a user-friendly computational tool that integrates clinical, biochemical, and immunological features, facilitating timely and accurate risk-stratification of patients with severe COVID-19. This will enable clinicians to make informed clinical decisions upon admission, predicting either recovery or clinical deterioration and promptly intervening when necessary. Independent, continuous validation of our computational tool in prospective settings and larger clinical trials will confirm its predictive accuracy, effectiveness, and usefulness in guiding therapeutic decisions.