TSCytoPred: a deep learning framework for inferring cytokine expression trajectories from irregular longitudinal gene expression data to enhance multi-omics analyses

Data collection and preprocessing

Cytokine and gene expression profiles, along with clinical information from COVID-19 patients, were obtained from the Clinical and Omics Data Archive (CODA) under accession number CODA_D23017 (Jo et al., 2022). This dataset contains longitudinal multi-omics data from 300 COVID-19 patients, including single-cell Ribonucleic acid (RNA) sequencing, cytokine profiling, and additional information such as the National Early Warning Score (NEWS) severity score. Additionally, it provides pseudo-bulk RNA sequencing data, which aggregates read counts across different cell types to generate patient-specific gene expression profiles. These aggregated profiles were used as the gene expression dataset in this study.

For analysis, we selected patients who had both cytokine and gene expression data available for three time points within a 15-day period (approximately two weeks). This timeframe was chosen based on reports that COVID-19 symptoms typically last up to two weeks (Raveendran, 2020; Wu, Yu & Lee, 2022). After filtering, a total of 93 patients met these criteria. Notably, time-series samples were collected at irregular intervals, meaning that the time gaps between successive time points varied across patients. The summaries of clinical severity, based on the distribution of NEWS scores for the 93 patients, and the distribution of time gaps between consecutive timepoints for each patient are provided in Supplementary Material S2.

Data preprocessing followed a similar methodology to previous studies analyzing cytokine profiles associated with COVID-19 disease progression (Han et al., 2024). Initially, only cytokines and genes shared across all samples were retained. Features with more than 50% missing values were removed to minimize bias, and the remaining missing values were imputed using a random forest imputation approach. To ensure consistency across samples, library size normalization was performed, followed by a log transformation to stabilize variance. After preprocessing, a total of 16,116 genes and 185 cytokines remained in the final dataset for analysis.

Selecting genes associated with cytokines

To prevent overfitting caused by high dimensionality while effectively capturing interactions between features for cytokine expression inference, the gene selection step follows a three-step process. First, genes that are known to be related to cytokines as transcription factors are selected based on the CytReg database (Carrasco Pro et al., 2018; Santoso et al., 2020). This database contains an extensive cytokine-gene network derived from functional assays, in vivo binding, and in vitro DNA binding, extending the existing literature on cytokine interactions. Next, TSCytoPred identifies additional genes associated with the previously selected cytokine-related genes based on protein-protein interactions from the STRING database (Szklarczyk et al., 2023). A combined score of 0.9 was used as the threshold for interaction significance. Finally, genes that are highly correlated with cytokines are selected. The Spearman correlation was calculated between each gene and its corresponding cytokine, and the top 50 genes with the highest correlation were retained. To enhance transparency, we further summarized the strength and significance of these correlations in Supplementary Material S3 (distribution of correlation coefficients across all cytokine–top 50 gene pairs) and Supplementary Material S4 (summary statistics for correlation and corresponding p-values including median, interquartile range, minimum, and maximum values for the top 50 genes per cytokine). As a result, our gene selection step identified 3,429 genes, which were then used as inputs for predicting time-series cytokine expression.

Inferring time-series cytokine expression based on deep learning

TSCytoPred is a deep learning model with an RNN-free architecture designed for multivariate time-series inference (Fig. 1). It is particularly suited for longitudinal datasets with a small number of time points and irregular time gaps. Let n represent the number of genes in the input and m represent the number of cytokines to be inferred. Let T = (t₁, t₂, …, t_k) denote an increasing sequence of times, where k is the total number of time points in the longitudinal data. The gene expression and cytokine expression observations at time t_i are represented as $x_i^g\in {\mathbb{R}}^n$ and $x_i^c\in {\mathbb{R}}^m$, respectively. The complete gene expression and cytokine expression matrices are denoted as $x^g=\left(x_1^g,\dots ,x_k^g\right)\in {\mathbb{R}}^{n\times k}$ and $x^c=\left(x_1^c,\dots ,x_k^c\right)\in {\mathbb{R}}^{m\times k}$, respectively. Additionally, the time gap vector, δ ∈ ℝ^k, is defined as the time lag between consecutive time points. Specifically, the time gap for the ith time point is given by δ_i = t_i − t_i−1.

Given the inputs x^g, x^c, and δ, TSCytoPred first extracts features using the MLP block l. This block consists of three fully connected (FC) layers with leaky ReLU activations, which produce the expansion coefficient predictors θ_ℓ as follows: (1)${\displaystyle h_{\ell ,1}=\text{leakyReLU}\left(W_{\ell ,1}{x}^g+b_{\ell ,1}\right),}$ (2)${\displaystyle h_{\ell ,2}=\text{leakyReLU}\left(W_{\ell ,2}{h}_{\ell ,1}+b_{\ell ,2}\right),}$ (3)${\displaystyle {\theta }_{\ell }=\text{leakyReLU}\left(W_{\ell ,3}{h}_{\ell ,2}+b_{\ell ,3}\right),}$ where W_ℓ,j for j ∈ 1, 2, 3 are the weights, and b_ℓ,j for j ∈ 1, 2, 3 are the bias terms for the MLP block ℓ. These expansion coefficients are then passed to the temporal interpolation block, which infers the cytokine expression for the given time points while accounting for the time gaps δ between them. In this framework, “time-series data” refers to longitudinal measurements from the same subject across multiple time points. For each subject, all available time points are processed jointly rather than independently, thereby preserving temporal dependencies across observations. The interpolation block then generates subject-specific cytokine trajectories over time, rather than a single static output, reflecting both the underlying gene expression dynamics and the irregular gaps between samples.

To explicitly address the irregular sampling intervals in our dataset, TSCytoPred incorporates δ directly into the temporal interpolation block. This design allows the model to use the actual elapsed time between successive samples to guide interpolation, rather than assuming uniformly spaced intervals. By doing so, patients with unevenly collected measurements are not excluded; instead, the model reconstructs continuous cytokine trajectories that respect the heterogeneous time gaps across individuals.

TSCytoPred employs linear interpolation to generate cytokine expression trajectories between time points t₁ and t_k, chosen for its computational efficiency and robust performance with sparsely sampled longitudinal data. Linear interpolation has demonstrated comparable or superior predictive performance in time-series tasks while requiring significantly less computational effort than alternatives such as nearest-neighbor methods (Lepot, Aubin & Clemens, 2017; Choi et al., 2023; Challu et al., 2023). To validate this design choice, we developed a TSCytoPred variant using cubic spline interpolation and conducted 5-fold cross-validation experiments. Linear interpolation consistently outperformed cubic spline interpolation across all evaluation metrics (Supplementary Material S5). While TSCytoPred can accommodate higher-degree methods such as polynomial, Hermite, or spline interpolation, these approaches typically require more than four time points, making them impractical for datasets with limited temporal sampling. The interpolation process generates continuous cytokine expression profiles across the entire time interval from t₁ to t_k, providing cytokine level estimates even at intermediate time points where gene expression data were unavailable (Fig. 1). The interpolation is computed as follows: (4)${\displaystyle {\hat{x}}_i^c={\theta }_{\ell ,1}+\frac{{\delta }_i}{t_k-t_i}\cdot \left({\theta }_{\ell ,k}-{\theta }_{\ell ,1}\right),}$ where ${\hat{x}}_i^c$ denotes the inferred cytokine expression at intermediate time point t_i, θ_ℓ,1 and θ_ℓ,k represent the inferred cytokine expression values at initial time t₁ and the final time point t_k, respectively, and δ_i represents the elapsed time gap from the initial time point t₁.

To train the model, the MAE loss was computed between the actual and inferred cytokine expression values for the non-missing time points. TSCytoPred was trained using the adaptive optimization algorithm Adam, with a learning rate of 10⁻⁵ with batch size of 30, and run for 5,000 epochs. Training was halted if the MAE loss did not improve over the course of 100 consecutive epochs. The model was implemented using the PyTorch (Version 1.6.0) (Paszke et al., 2019), Scikit-learn (Version 1.7.1) (Pedregosa et al., 2011), Pandas (Version 2.3.1) (Wes McKinney, 2010), and Numpy (Version 2.2.6) (Harris et al., 2020). Hyperparameters were tuned using grid search, with each experiment repeated five times. The hyperparameters that minimized the average MAE loss were selected as optimal values. The hyperparameter optimization results including the number of hidden nodes, number of layers, learning rate, and batch size are summarized in Supplementary Material S6. The MLP block consisted of three fully connected (FC) layers, with 1024, 512, and a final layer size equal to the number of cytokines to be inferred. The model was evaluated on a single-node server equipped with 40 Intel(R) Xeon(R) Silver 4114 CPU cores running at 2.20 GHz, 250 GB of main memory, and an NVIDIA RTX 8000 GPU with 48 GB of memory.

Performance evaluation of TSCytoPred for inferring time-series cytokine expression

TSCytoPred was benchmarked against several widely used regression models, including ML approaches—ElasticNet, Lasso, Ridge, and Linear Regression—as well as the statistical model ARIMA. We further compared TSCytoPred with a neural network-based regression model (NN), which treated each time point as a separate sample and used a similar MLP block-based architecture with the same number of nodes and fully connected layers as in our model. In addition, we evaluated recurrent neural network models, including a LSTM network and a hybrid convolutional neural network using 1D-CNN combined with LSTM (CNN-LSTM), to assess their inference performance. The ML regressors were implemented using ‘Scikit-learn’ python package. A 5-fold cross-validation was performed at the patient level (i.e., each patient’s three time points were kept within the same fold to prevent any data leakage across folds), and the hyperparameters were optimized by grid search based on the average MAE (see Supplementary Material S7). The optimized hyperparameter settings for each regression model were as follows: ElasticNet (α = 0.1), Ridge (α = 10), Lasso (α = 0.1); LSTM (two layers with 2,000 and 1,000 hidden units, respectively); and CNN-LSTM (kernel size = 3, padding = 1, 2,000 channels, followed by two LSTM layers with 1,000 and 250 hidden units). To ensure a fair comparison, the same input data—comprising the genes selected by our gene selection step—was used for all comparison methods. Model performance was evaluated using multiple metrics. MAE and root mean square error (RMSE) quantify the magnitude of prediction errors, with RMSE placing greater weight on larger errors. The R² coefficient indicates the proportion of variance in cytokine levels explained by the model. Mean absolute percentage error (MAPE) expresses prediction errors as a percentage of actual values, enabling scale-independent comparisons across cytokines. Spearman correlation (CORR) measures the rank-based association between predicted and actual cytokine levels, capturing the consistency of the prediction ordering even when absolute values differ. Evaluations were conducted using multi-omics COVID-19 datasets containing time-series gene expression and cytokine profiles collected from 93 patients across three time points within a 15-day period. A detailed description of this dataset can be found in ‘Methods’ section.

As summarized in Table 1 (with 95% confidence intervals provided in Supplementary Material S8), TSCytoPred outperformed the other models in cytokine expression inference, achieving the highest average R² coefficient of 0.257, along with the lowest average MAE (0.437) and MAPE (0.118). ElasticNet achieved the second-best R² of 0.247, while the NN model had the second-best performance in the other metrics, with an MAE of 0.442 and the best MAPE of 0.116. By contrast, both LSTM and ARIMA produced negative R² values and higher MAE, indicating that, for datasets with limited and irregularly sampled time points, recurrent models and statistical approaches are disadvantaged due to their reliance on regularly spaced temporal measurements and their need for longer sequences to achieve stable training. Given these poor performances, we excluded LSTM and ARIMA from further experiments. To further evaluate robustness, we also conducted a 3-fold cross-validation, which yielded results consistent with the 5-fold setting, confirming that TSCytoPred’s performance remains stable under different partitioning schemes (Supplementary Material S9). Furthermore, TSCytoPred demonstrated greater stability in performance across the 5-fold cross-validation, exhibiting less variation compared to the other methods. We recognize, however, that the overall R² values remain relatively modest, suggesting limited variance explained. To investigate this further, we conducted sensitivity analyses by removing cytokines with varying proportions of outlier samples (defined using interquartile range (IQR)-based thresholds of >5%, >3%, and >1%). These additional results (Table 2, Supplementary Material S10) showed that after removing outlier cytokines, the average R² of TSCytoPred improved from 0.257 to 0.362, while most comparison methods also achieved R² values above 0.3. Importantly, TSCytoPred remained the best-performing method across all scenarios, supporting the robustness of the model even when noisy cytokines were excluded.

To further assess robustness, we also repeated the evaluation using a larger subset of COVID-19 patients with two time points within 9 days (125 patients), and TSCytoPred continued to outperform other models in most metrics (Supplementary Material S11). In addition, we tested the model on the original three–time point cohort by predicting cytokine abundance dynamics in reverse, where TSCytoPred again outperformed baseline methods, including NN and CNN–LSTM (Supplementary Material S12). These results confirm that TSCytoPred provides reliable inference performance even when longer-term cytokine trajectories are unavailable or when different temporal dynamics are considered.

We also assessed model performance by comparing actual and inferred cytokine expression values through visualization. For this purpose, the dataset was randomly split into an 8:2 ratio for training and testing solely for visual evaluation, with models trained on the training set. This split was not used for reported performance metrics; instead, all quantitative evaluations in this study were conducted using five-fold cross-validation to ensure robustness and reproducibility. We first visualized the cytokine expression profiles using uniform manifold approximation and projection (UMAP) at the sample level to assess how closely regression models infer to actual data, categorized by clinical COVID-19 risk (Fig. 2). Ideally, inferred profiles (triangles) should cluster closely with actual measured profiles (circles) within each clinical risk category. TScytoPred exhibited the strongest alignment, as indicated by the close proximity of inferred to actual data points, especially notable within the medium-risk category (green). In contrast, methods such as Lasso displayed significant discrepancies, with inferred points broadly dispersed from actual measurements. Similarly, Ridge and Linear methods showed noticeable separations between inferred and actual profiles, reflecting lower accuracy in inference.

To further evaluate the performance, we examined the individual cytokine predictions. The predicted cytokine expression values in the test data were compared to the actual values using the Wilcoxon signed-rank test. All p-values from these tests were adjusted for multiple comparisons using the false discovery rate (FDR) correction (Benjamini–Hochberg procedure). Among the 185 cytokines, 84.32% (156 cytokines) of those inferred by TSCytoPred showed no significant difference from the actual values (p-value > 0.05) as shown in Table 3. In contrast, other regression methods exhibited significant discrepancies between predicted and actual cytokine expression, with the following numbers of cytokines showing no significant differences: 152 (NN), 149 (ElasticNet), 155 (Linear), 147 (Lasso), 153 (Ridge) and and 80 (CNN-LSTM). As an example, Supplementary Material S13 presents boxplots comparing the distributions of actual and inferred cytokine expression for a few key cytokines associated with COVID-19. TSCytoPred exhibited the closest alignment to the actual expression distribution, effectively capturing the patterns observed in the real data. This demonstrates that TSCytoPred outperforms the other methods in inferring individual cytokine patterns.

Effectiveness of gene selection in TSCytoPred

TSCytoPred employs a comprehensive gene selection strategy incorporating three key steps to identify gene sets associated with cytokines, facilitating accurate inference of cytokine expression. To assess the effectiveness of each step, we developed three variants of TSCytoPred by sequentially modifying or removing components of the gene selection procedure. Specifically, the ‘CytReg’ variant selects only transcription factor genes directly associated with cytokines based on the CytReg database. The ‘CytReg + STRING’ variant further includes additional genes connected to these cytokines through protein-protein interactions from the STRING database. Finally, the complete approach, labeled ‘CytReg + STRING + CORR’, represents the full TSCytoPred method, where we evaluated performance using varying numbers of top-ranked genes according to Spearman correlation with each cytokine (e.g., ‘top 10’ refers to selecting the 10 genes with the highest correlation per cytokine). Additionally, we tested another variant utilizing mutual information (MI)—a common alternative feature selection method in time-series analysis—in place of Spearman correlation. Performance across all variants was evaluated using 5-fold cross-validation.

The comparative results, summarized in Table 4 (with 95% confidence intervals provided in Supplementary Material S14), demonstrate that the full TSCytoPred approach, integrating genes selected from CytReg, STRING, and Spearman correlation rankings, achieved the highest predictive accuracy, with an average R² of 0.257 and MAE of 0.437. In contrast, replacing correlation-based selection with mutual information did not yield significant improvements, despite using a comparable number of genes. These findings underline the effectiveness and importance of our gene selection strategy for enhancing cytokine expression inference from gene expression data.

Improvement of the severity prediction through cytokine expression inference using TSCytoPred

Cytokine levels have been widely recognized as predictors of mortality and illness severity in COVID-19 infections (Del Valle et al., 2020; Cabaro et al., 2021; Onuk et al., 2023). In this study, we explored whether utilizing cytokine data inferred through TSCytoPred could enhance severity outcome prediction, specifically using the NEWS. The NEWS system assesses the severity of illness based on seven criteria: respiratory rate, oxygen saturation, temperature, blood pressure, and heart rate (Martín-Rodríguez et al., 2022). The score ranges from 0 to 21, with three points awarded per criterion, where higher scores indicate more abnormal values and a more severe condition. In our COVID-19 dataset, clinical information for the NEWS severity score was available for 113 patients. Among them, 93 had both gene expression and cytokine data at three time points, while the remaining 20 had only gene expression data. A detailed description of the dataset can be found in the ‘Methods’ section.

To investigate whether inferred cytokine data could enhance the prediction of disease severity outcomes, we compared regression models trained with and without inferred cytokine profiles. We employed an ElasticNet regressor—selected based on its robust performance in previous evaluations—and trained it using available cytokine data along with inferred cytokine expression values for 20 patients who originally had only gene expression data. Cytokine data inferred by TSCytoPred and other baseline methods were evaluated and compared against predictions from a regressor trained solely on existing cytokine measurements (i.e., without inferred cytokine data).

The average R² and MAE from 5-fold cross-validation showed that training the regressor with the inferred cytokine profiles from TSCytoPred led to a 6.12% performance improvement, with the average R² increasing from 0.294 to 0.312 (Fig. 3). In contrast, other methods showed little improvement, with the NN model achieving the second-best R² of 0.297. These findings suggest that TSCytoPred could help the regression model to effectively leverages inferred cytokine data to enhance the prediction of clinical severity outcomes, especially for patients lacking direct cytokine measurements.

Enhancing clinical risk prediction for COVID-19 severity outcome using cytokine profiles

While predicting severity scores for each timepoint is important, directly predicting a patient’s overall clinical risk based on cytokine profiles could offer a more efficient approach for practical applications. Ultimately, patients will be classified into clinical risk categories based on the NEWS score threshold, and their responses will be predicted accordingly. In this experiment, we transformed the NEWS severity score values into categorical outcomes based on the predefined thresholds used in hospitals (Welch, Dean & Hartin, 2022), as shown in Table 5. A support vector machine (SVM) classifier was selected and constructed using the ‘Scikit-learn’ package with default parameters due to its superior classification performance compared to methods including decision trees and Naive Bayes (Zhang et al., 2017), taking cytokine expression profiles as input.

The performance of the classifier for clinical risk prediction was assessed by training the models with cytokine data inferred by TSCytoPred for the 20 patients lacking direct cytokine measurements, along with data inferred by other baseline methods. Model performance was evaluated using F1-scores based on 5-fold cross-validation, under two distinct prediction scenarios. In the first scenario (time-point-level prediction), models aimed to predict clinical risk at individual time points. In the second scenario (patient-worst-level prediction), the goal was to predict the highest clinical risk level that a patient might experience throughout the entire course of illness. Predicting the worst-case clinical risk is especially valuable for clinicians, as it allows them to anticipate patient needs and implement targeted interventions. For example, if a patient exhibited clinical risk levels of ‘Low,’ ‘Medium,’ and ‘Low’ across three different time points, the task would be to correctly identify ’Medium’ as the highest risk level experienced by that patient.

As summarized in Table 6 (with 95% confidence intervals provided in Supplementary Material S15), the addition of inferred cytokine samples from TSCytoPred improved the average F1-score in both scenarios. For timepoint-level prediction, the average F1-score improved from 0.565 to 0.602, while for patient-level prediction, it improved from 0.568 to 0.592. TSCytoPred outperformed all other methods in both severity score and clinical risk predictions. These results demonstrate that TSCytoPred is an effective tool for inferring missing cytokine profiles that would otherwise be discarded, and that the inclusion of inferred data significantly enhances model performance.