{'fraction_non_alphanumeric': 0.03994403499783234, 'fraction_numerical': 0.02256335474717219, 'mean_word_length': 5.24271128060032, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-omics personalized network analyses highlight progressive disruption of central metabolism associated with COVID-19 severity

The clinical outcome and disease severity in coronavirus disease 2019 (COVID-19) are heterogeneous, and the progression or fatality of the disease cannot be explained by a single factor like age or comorbidities. In this study, we used system-wide network-based system biology analysis using whole blood RNA sequencing, immunophenotyping by flow cytometry, plasma metabolomics, and single-cell-type metabolomics of monocytes to identify the potential determinants of COVID-19 severity at personalized and group levels. Digital cell quantification and immunophenotyping of the mononuclear phagocytes indicated a substantial role in coordinating the immune cells that mediate COVID-19 severity. Stratum-specific and personalized genome-scale metabolic modeling indicated monocarboxylate transporter family genes (e.g., SLC16A6), nucleoside transporter genes (e.g., SLC29A1), and metabolites such as α-ketoglutarate, succinate, malate, and butyrate could play a crucial role in COVID-19 severity. Metabolic perturbations targeting the central metabolic pathway (TCA cycle) can be an alternate treatment strategy in severe COVID-19.

Correspondence ujjwal.neogi@ki.se In brief Ambikan et al. used blood cell transcriptomics, immunophenotyping, plasma metabolomics, and single-celltype metabolomics of monocytes to identify the system-level metabolic rewiring in COVID-19 patients. Integrative omics improved the clinical definition of the risk group of COVID-19 severity. The personalized and group-specific metabolic models indicated the essential role of transporters and metabolites of central metabolism (TCA cycle) in COVID-19 severity. This can lead to an alternate treatment strategy through metabolic perturbations of central metabolism in severe COVID-19.

INTRODUCTION
Coronavirus disease 19 , caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), presents mild to moderate symptoms in most of the infected patients. Several studies, including ours, have indicated altered systems-level metabolic disruption in SARS-CoV-2-infected patients compared between metabolic and signaling molecular effectors driving disease severity has not been systematically examined. Further, in addition to comorbidities, the high physiological heterogeneity has not been considered when addressing the mechanisms of population-level predisposition to the severe phenotype.
Viruses are obligatory parasites and entirely rely on their hosts for replication. This reliance is evidenced by experimental findings showing significant metabolic flux disruption in the host cells. This is further highlighted by observed variations in the cell-specific viral replication and production, leading to changes in the host metabolism (Yu et al., 2011). Virus infection leads to significant metabolic alterations in the host, such as altered glycolysis rates and changes in ATP production rates (Aller et al., 2018). Therefore, the changes in energy metabolism can be seen as an evolving property of the combined host-virus metabolic system. This could be related to changes in host cellular demands arising from viral production, increased demand from immune cell activation, and damage due to toxicities (e.g., oxidative damage) (Molenaar et al., 2009). Even though the primary site of infection is the upper respiratory tract, SARS-CoV-2 can invade several organs, tissues, and cells of the body (Delorey et al., 2021). Therefore, a systems-level characterization of metabolic and signaling dysregulation can comprehensively and quantitatively capture the implications of COVID-19 infection throughout the body.
Systems biology approaches have shown promise in systematically elucidating complex disease-underlying mechanisms by integrating different layers of omics information. The application of the methods, including machine learning, biological network analyses, and multi-omics data integration, has allowed us to generate multi-omics integrated networks to characterize the host response to SARS-CoV-2 . Interrogating these integrated networks permits identifying and understanding the sequence of specific biological elements coordinately altered in pathological conditions. Further, comprehensive metabolic characterizations allow us to simulate system-level metabolic fluxes through genome-wide metabolic modeling in human cells and tissues . Combining biological network analysis and genome-scale metabolic modeling (GSMM) as multiomics integrative approaches can attain holistic and temporally dynamic characterizations of complex rearrangement in response to viral infection and can computationally predict potential therapeutic targets. Importantly, these characterizations permit a group-specific and personalized characterization of global metabolism in healthy and progressive disease states.
In the present study, we employ data-driven multi-omics, whole blood transcriptomics, and plasma metabolomics (Krishnan et al., 2021) characterization and clustering of patient data through the construction of integrated networks. Further, the multidimensional stratification of patients through similarity network fusion (SNF) enables addressing population heterogeneity and redefining the clinical classification of COVID-19 severity. We also develop a GSMM encompassing both host and viral metabolic demands that permit metabolic characterization at a personalized level and attain a data-driven patient stratification. Subsequently, a combination of metabolic flux prediction, topological network analysis, reporter metabolite analysis, and in silico target prediction of the personalized networks permits identifying commonalities and patient-specific targets and hubs. Finally, we measure the intracellular metabolites that are predicted to be the key regulators during SARS-CoV-2 infection in the specific cell type identified by digital cell quantification (DCQ), immunophenotyping, and single-cell RNA sequencing (scRNA-seq) data (obtained from Zhang et al., 2020) from COVID-19 patients. We combine integrated network analysis and personalized GSMM for COVID-19 patients at the metabo-transcriptomic levels that could pave the way for personalized treatment strategies dependent on the patient's metabolic state.

Study design and patient cohorts
Based on the oxygen (O 2 ) requirement and the clinician's decision, we clinically categorized COVID-19 patients into hospitalized-mild (n = 26, O 2 requirement < 4 L/min) and hospitalized-severe (n = 11, O 2 requirement > 4 L/min). The patient characteristics were presented elsewhere (Krishnan et al., 2021). Additional COVID-19 PCR-negative samples were also collected (healthy control [HC], n = 31) and tested for SARS-CoV-2 antibody levels. The serology identified 10 participants as SARS-CoV-2 antibody-positive whose samples were termed as convalescent. RNA sequencing (RNA-seq) was performed on the whole blood samples obtained from SARS-CoV-2-infected and -noninfected patients. We used DCQ to quantify the different cell types, followed by differential analysis between clinically identified groups after adjusting for the high abundant cell types. Immunophenotyping using flow cytometry further validated the changes in mononuclear phagocytes (MNPs). To improve the classification of clinical COVID-19 severity, transcriptomics and metabolomics data were integrated using SNF, leading to the re-classification of the patient groups into healthy, mild/moderate, and severe. SNF-classified patient groups led to contextspecific GSMM both at personalized and group levels. The flux balance analysis (FBA) was performed to study the metabolic rearrangement in deeper resolution and identify the essential genes and metabolites associated with SNF-defined COVID-19 severity in groups. Previously published scRNA-seq data  was re-analyzed to examine the expression of the essential genes that are associated with disease severity in specific cell types. Finally, single-cell-type metabolomics analysis (sctMetabolomics), targeting TCA-cycle intermediates, confirmed the metabolic rearrangement in the monocytes. The complete study design is presented in Figure 1.
Distinct changes in immune cell types in hospitalizedsevere patients The transcriptomics data was generated from whole blood samples; therefore, the gene expression variation may be due to altered gene regulation in cell types whose abundances differed between patient groups. Thus, we have performed DCQ that computationally estimates the cell type proportions from bulk blood RNA-seq data. DCQ uses a cell-specific expression profile of 18 blood cell types obtained from the Human Protein Atlas (Uhlen et al., 2017) and a deconvolution algorithm adapted from the method estimation of proportions of immune and cancer cells (EPIC) (Racle and Gfeller, 2020). Proportions of 18 blood cell populations in all the samples were computed and analyzed for groupwise differences. Kruskal-Wallis's test identified four cell types, namely classical monocytes (CM), neutrophils, mucosa-associated invariant T (MAIT) cells, and gd-T cells, which differed among the patient groups (adj. p < 0.05; Figure 2A). The pairwise analysis identified that in both hospitalized-mild and -severe groups, neutrophils and CM displayed significantly increased abundance in COVID-19 patients compared with HC, whereas gd-T cells, memory B cells, and MAIT cells had lower abundance in severe COVID-19 patients compared with HC (adj. p < 0.05; Figure 2B). The pairwise co-expression landscape of marker genes of specific cell types showed a preserved correlation between the signature genes in the immune cell types in HC ( Figure 2C) but a dysregulated correlation in convalescent ( Figure 2D), mild ( Figure 2E), and severe ( Figure 2F) patients. Specifically, the associations between marker genes of T cells and dendritic cells (DC) and the associations between granulocytes (neutrophils/eosinophils) and DC or T cells were altered, indicating an impaired interaction between these cell types during COVID-19 infection at a system level ( Figure S1A) that may lead to an inability of the innate immune cells to elicit an effective adaptive immune response in severe COVID-19 patients. The GSEA performed on the lost marker genes associations indi-cated significant changes in PI3K-Akt signaling, chemokine signaling pathway, NF-kB, and MAPK signaling pathways ( Figure S1B).

Immunophenotyping of blood MNPs associated with COVID-19 severity
The loss in the co-expression of the marker genes between the T cells and myeloid cells observed in the transcriptomics data prompted us to characterize these populations in peripheral blood mononuclear cells (PBMCs) of COVID-19 mild and severe patients. We performed immunophenotyping of MNPs consisting of monocytes and DCs and innate regulatory cells such as myeloid-derived suppressor cells (MDSCs) and lowdensity granulocytes (LDG). We used CD4 + and CD8 + T cells from our earlier study (Krishnan et al., 2021). The gating strategy and markers used to define the different cell populations of MNPs are provided ( Figure S2; Table S1), and representative contour plots for HC, mild, and severe patients are shown in Figure 3A. Integrating all analyzed markers and projection of the gated populations into the UMAP space revealed differences in several cell populations between HC and COVID-19  Figures 3C and  S3B. A decrease in G-MDSC and an increase in M-MDSC (albeit very low levels) and LDG frequencies were observed in COVID-19 patients, with differences being more pronounced in patients with severe disease ( Figure 3C). In the DC and monocyte populations, significant decreases in relative frequencies between HC and COVID-19 patients were noted for plasmacytoid DC (pDC), CD16 + non-classical monocytes, which also includes monocyte-derived DCs (NCM/MD-DC), and total CD16 À classical DC (cDC). On the contrary, the frequencies of the CD303 À DC5 population (Villani et al., 2017) and of the CM were increased in COVID-19 patients compared with HC. We also observed a significant decrease in CD16 + CD141 low CD1c À DC4 population in mild COVID-19 patients compared with HC ( Figure 3C). Although the total CD16 À cDCs were decreased, no difference in the relative frequency of CD141 + DC1 or CD1c + DC2/DC3 populations (Figure S3B) was observed. In line with the DCQ, the association between the T cells (CD4 + and CD8 + ) and some of the MNPs (e.g., pDC) was altered ( Figure 3D). We then assessed the expression of the chemokine (C-C motif) receptor type 2 (CCR2), CX3 chemokine receptor 1 (CX3CR1), and C-C chemokine receptor type 5 (CCR5) in the MNPs, using flow cytometry ( Figures 3E-3G). Our analysis showed an increase in the frequency of CCR5 + CM, IM, NCM/MD-DC, DC4/NCM cells, and DC1 cells, especially in the mild COVID-19 patients (Figure 3G). The frequency of CCR5 + G-MDSCs was significantly decreased in COVID-19 patients in a severity-dependent way. Similarly, CCR2 + IM, NCM/MD-DC, DC4, and DC5 were increased in frequency, and CCR2 + DC1 cells were decreased in COVID-19 patients, compared with HC ( Figure 3G). We also observed a significant decrease in the frequency of CX3CR1 + pDC (both in mild and severe , and in CM (in severe COVID-19 patients) and NCM/MD-DC (in mild COVID-19 patients), compared with HC ( Figures 3G and S3C). The median fluorescence intensity (MFI) of the chemokine receptors was also altered in COVID-19 patients, with altered CCR5 expression more predominant in mild patients and altered CCR2 expression in all COVID-19 patients, compared with HC ( Figures 3E and S3D). These data highlight essential changes in the MNP landscape in the blood of COVID-19 patients, including a phenotype consistent with chemokine-mediated migration of MNPs (particularly of DCs) in the infected lung and additionally defective interaction of MNP with lymphocytes, which probably prevent efficient viral clearance.
System-level whole blood transcriptomics signature differentiates COVID-19 patients To identify the system-level host response following the SARS-CoV-2 infection, we performed a transcriptomic profile of the total RNA isolated from whole blood by RNA-seq. A UMAP projection of the data showed a clear separation of the COVID-19 patients from HC and convalescent controls but failed to efficiently distinguish between hospitalized-mild and -severe patients ( Figure 4A). Differential gene expression (DGE) analysis (adj. p < 0.05 and log 2 -fold change R1.5) identified 581 genes that were differentially regulated between the HC and COVID-19 patients ( Figure 4B; Data S1) and 154 genes differentially regulated between hospitalized-mild and -severe patients (Figure 4C; Data S1). While comparing the HC with the hospitalized-mild and -severe groups, 445 and 1,068 genes were significantly different between the two groups, respectively (Data S1). The hierarchical clustering analysis (HCA) based on the significantly regulated genes also identified distinct patient clustering but with a heterogeneity among the COVID-19 patients (Figure 4D). The gene set analysis between the HC and the COVID19 patients identified distinct upregulations of pathways related to antiviral response (e.g., NOD-like receptor signaling, RIG-I-like signaling, and TNF signaling pathways), metabolism (e.g., fructose mannose metabolism, oxidative phosphorylation (OXPHOS), and pentose phosphate pathway), and pathways like C-type lectin receptor signaling pathway, complement, and coagulation cascades that are related to thromboembolism ( Figure 4E). Interestingly, a comparison between mild and severe patients' profiles showed that pathways like PI3K-Akt, mTOR, AMPK, and HIF-1 signaling were distinctly upregulated in severe patients ( Figure 4E). These pathways are the master regulators of central carbon and energy metabolism (e.g., TCA cycle, OXPHOS, and pyruvate metabolism) and are regulated by the SARS-CoV-2 replication Codo et al., 2020;Krishnan et al., 2021). The gene set analysis result for group-specific (mild and severe) comparisons with HC is provided ( Figures 6 and S4A).
Data-driven patient stratification results in the reclassification of COVID-19 subgroups Clinical categorization of COVID-19 patients relies on their oxygen requirement. However, our previous study identified that O 2 need at hospitalization did not predict mortality . Moreover, the HCA analysis identified heterogeneity among the significantly dysregulated genes in COVID-19 patients ( Figure 4D). Therefore, we sought to improve the classification of clinical COVID-19 severity definitions as driven from a molecular point of view through a data-driven integrative (C-F) Network visualization of co-expression among marker genes of various cell types in samples of healthy control (HC) cohort (C), convalescent cohort (D), hospitalized-mild cohort (E), and hospitalized-severe cohort (F). Each node in the network represents marker genes of the corresponding cell type, and node size is relative to the mean expression value (TPM). Edge denotes a significant Spearman correlation (adj. p < 0.001) between marker genes in HC, convalescent, and hospitalized-mild cohorts while adjusted p < 0.2 in hospitalized-severe.  Article approach capable of handling the high patient-patient molecular heterogeneity. We integrated information from both transcriptomics and metabolomics through SNF, which has been shown to capture both shared and complementary information from different data sources and produce clinically relevant subtypes in multiple studies . Through SNF, we observed that while individual omics did not show evident similarity patterns ( Figures 5A and 5B), the integrative characterization of the samples pointed to four clusters ( Figure 5C; Table S2). Among these, we observed that many HC (52.4%, 11/21;) samples shared high similarities with convalescent (60%, 6/10) samples (SNF cluster 4; SNF-4) ( Figures 5C-5E), in line with our observations at transcriptomic-only level ( Figure 4A). Furthermore, 50% of hospitalized-severe samples shared high similarity with hospitalized-mild samples (SNF-3). In addition, we observed that SNF-2 included samples from all four groups, including HC (38%, 8/21), convalescent (40%, 4/10), hospitalized-mild (6.9%, 2/29), and hospitalized-severe (27%, 3/11) samples, suggesting that homogeneous features exist in each layer of data type. Furthermore, SNF-1 was dominated by hospitalized-mild samples (69.0%, 20/29) ( Figures 5C and 5D). Based on the plasma inflammation profile identified by the interleukins (IL) IL6, IL7, IL8, IL10, IL15, and IL18 ( Figure 5F)-measured using the Olink Immuno-oncology panel (Krishnan et al., 2021)-we defined SNF-1 as COVID-19 mild/moderate, SNF-3 as COVID-19 severe, and SNF-4 as HC. In subsequent analyses, we used these newly defined clusters except for SNF-2 due to its mixed populations that may share common molecular features.
Multi-omics network topology analysis identifies severity-specific mechanisms Integrated networks based on transcriptomics and metabolomics data were generated to identify severity-specific omics hallmarks and associated mechanisms. Topological network analysis (adj. p < 10 À5 , Spearman rank correlation) among features from samples of SNF-1, 3, and 4 was used to create a network, followed by a community identification . The analysis highlighted nine communities (c1-c9) of co-expressed genes and metabolites ( Figure 5G). The expression direction of genes and metabolites in SNF-3 (severe) compared with SNF-1 (mild/moderate) was mapped to the network to define severity associations at the feature community level. Out of the total genes in c1, 42% of them were found significantly (adj. p < 0.05) upregulated in SNF-3 (severe) compared with SNF-1 (mild/moderate), and thus it was termed a severityspecific community. Among the rest of the genes in c1, 6.6% were upregulated in SNF-1. We employed centrality analysis to identify and rank the importance of genes and metabolites based on how connected they were to the network (i.e., degree centrality). We observed (Data S2) that c1 was the most central commu-nity, thus showing that the most connected molecular elements tended to be associated with coordinated expression changes related to COVID-19 severity state. Among the top 5% central elements in the network (542 nodes), 290 genes/metabolites were found in c1 and included the glycolytic genes ADP-dependent glucokinase (ADPGK), phosphoglycerate mutase 1 (PGAM1), and the TCA-cycle and OXPHOS gene succinate dehydrogenase complex iron-sulfur subunit B (SDHB). We observed that mannose, a biomarker for COVID-19 (Krishnan et al., 2021), was also found in c1. Enrichment analysis showed that c1 was associated with NOD-like receptor signaling, chemokine signaling, Fc gamma R-mediated phagocytosis and platelet activation pathways, and other pathways such as NF-k B, Notch, RIG-I, HIF-1, and FoxO signaling ( Figure 5H). The enrichment analysis of the other communities (c2-c9) is given ( Figure S5). These observations indicate that substantial signaling alterations in the immune response are closely linked with central metabolism through mannose, ADPGK, PGAM1, and SHDB, all of which are associated with disease severity.

Metabolic modeling highlights host metabolic dysregulation in COVID-19
Since prior omics analysis has identified metabolic and signaling pathway dysregulation in COVID-19 ( Figure 4E), we have used context-specific GSMM to find how the metabolic reactions are impacted upon COVID-19 and their cellular localization further to increase the resolution of perception on metabolic rearrangement. We have generated group-specific GSMMs based on the newly defined SNF clusters by integrating transcriptomics data and constrained the nutrient uptake of the models based on metabolomics data. Each model accounted for the intracellular energy requirements for the growth of host cells (see STAR Methods). Models were deemed feasible if they passed 57 biological functions . After FBA and based on metabolic flux distributions, we identified 100 reactions differentially active (with or without flux) or with opposite reaction directions among SNF clusters (Data S3). Further, we identified 15 transporter reactions with unique flux in SNF-1 (mild/moderate) and SNF-3 (severe) compared with SNF-4 (HC) ( Table S3). Among the transporter reactions, there were reactions involved in the transport of TCA-cycle intermediates such as cis-aconitate, a-ketoglutarate, succinate, malate, and fumarate between cytosol and mitochondria. Out of these 15 reactions, 13 had flux only in SNF-1 (mild/moderate) and SNF-3 (severe), indicating active transport of the metabolites in COVID-19 patients, which then potentially feed the TCA-cycle pathway ( Figure 6A). We used previously published single-cell transcriptomics data , which characterized patient groups similar to ours. The members of the mitochondrial carrier family (SLC25) (SLC25A1 and SLC25A11) were highly expressed   (F) Mean-scaled interleukins (IL) of plasma inflammation profile identified by the IL6, IL7, IL8, IL10, IL15, and IL18 to define the severity of the data-driven clusters. (G) Community characterization results obtained from topology analysis of the association network between genes and metabolites among samples of SNF-1, 3, and 4. Node size is relative to the centrality (degree) of the feature (gene/metabolite). Node colors denote log 2 -fold change values of the feature in SNF-3 (severe, n = 10) compared with SNF-1 (mild/moderate, n = 22). Each edge denotes a significant Spearman correlation (adj. p < 0.00001).
(H) Significantly enriched pathways (adj. p < 0.01) in the most central community found in (F) (community 1) based on average node centrality (degree).
(log expression > 2) in more cells in the monocyte populations ( Figures 6B and S6). As we observed that signatures of monocytes were enriched both in the DCQ (Figure 2) and in immunophenotyping analyses (Figure 3) and metabolism in monocytes plays an important role in regulating innate immune function (Pence and Yarbro, 2019), we performed sctMetabolomics targeting TCA-cycle intermediates in monocytes. We observed increased cellular levels of a-ketoglutarate, malate, and citrate in monocytes in COVID-19 samples ( Figure 6C). To further handle the individual heterogeneity underlying the molecular profiles of all patients, we generated personalized GSMMs to characterize COVID-19-specific and severity-specific metabolic alterations. FBA identified 274 metabolic reactions with distinct fluxes across patients ( Figure 6D; Data S3), many of which recapitulated the above-described group-specific observations. Reactions that had flux in the same direction in more than or half of the samples in SNF-1 and SNF-3 and opposite or no flux in more than or half of the samples in SNF-4 were regarded as COVID-19-specific reactions. By contrast, reactions that had flux in the same direction in more than or half of the samples in SNF-3 and opposite or no flux in more than or half of the samples in SNF-1 were regarded as COVID-19 severity-specific reactions. The analysis provided 16 reactions as COVID-19 specific and 10 reactions as COVID-19 severity specific ( Figure 6D; Data S3).

FBA-derived network topology analysis identified essential genes and metabolites in COVID-19
The results above indicated substantial rearrangement of intracellular metabolic fluxes concomitant with COVID-19 and its severity. To further explore the biological implications of the metabo-transcriptomic alterations while still considering the biologically driven characterization provided by the predicted flux dis-tributions, we performed a topological analysis of the personalized and group-specific metabolic networks. These weighted networks were created using reactions that displayed divergent metabolic fluxes, where nodes were reactants, products, and the associated enzymatic genes (Data S3), and edge weight was scaled down measurement of the reaction fluxes. Further, network communities were identified, and centrality was calculated to prioritize nodes (Data S3). We observed that 7 genes and 51 metabolites showed very high centrality (Figure 7A; betweenness > mean across samples; 83.9 for genes and 22 for metabolites). These included the monocarboxylate transporter SLC16A6 and nucleoside transporter SLC29A1 and metabolites such as a-ketoglutarate in the cytosol, succinate and malate in mitochondria and cytosol, and butyrate in the cytosol and extracellular space, which were identified as uniquely essential for COVID-19 samples. Topology analysis of SNF-3 (severe) networks also identified a-ketoglutarate as the most central among other intermediate products of the TCA cycle such as cis-aconitate, malate, fumarate, succinate, and mitochondrial transporter genes like SLC25A1 and SLC25A11 ( Figure 7B). The reactions associated with SLC25A1 and SLC25A11 displayed fluxes ( Figure 6D) only in SNF-1 (mild/moderate) and SNF-3 (severe) but not in any other SNF cluster, further reinforcing the rearrangement of malate and a-ketoglutarate metabolism observed at sctMetabolomic level ( Figures 6A  and 6B). Moreover, gene essentiality analysis identified mitochondrial genes as essential in COVID-19 patients ( Figure S7), including NADH: ubiquinone oxidoreductase core subunit S2 (NDUFS2) and NADH: ubiquinone oxidoreductase subunit B3 (NDUFB3) in OXPHOS, glycine N-methyltransferase (GNMT) in glycine metabolism, 5-oxoprolinase, ATP-hydrolysing (OPLAH) in glutathione metabolism, and amidohydrolase domain containing 1 (AMDHD1) in histidine metabolism, all of which were part of the central community c1 ( Figure 5G). To find the cell types where the top central genes are mostly expressed, we used the published scRNA-seq data . COVID-19specific essential gene SLC16A6 was mostly expressed in monocytes and monocyte-derived DCs ( Figures 7C and S8), whereas SLC29A1 was mostly expressed in NCM and proT cells ( Figures 7C and S8).
DISCUSSION This study combined system-level transcriptomics, DCQ, blood cell immunophenotyping of MNPs, integrated network analysis, and generated personalized and group-specific GSMM of COVID-19 patients to provide comprehensive systems-level characterizations in SARS-CoV-2 infection. The transcriptomics data identified distinct upregulation of pathways related to the antiviral response and metabolism in the COVID-19 patients. While comparing the mild and severe COVID-19 patients', pathways that are master regulators of central carbon and energy metabolism (e.g., TCA cycle, OXPHOS, pyruvate metabolism, etc.), the PI3K-Akt, mTOR, AMPK, and HIF-1 signaling were distinctly upregulated, which are also known to be regulated during SARS-CoV-2 infection Krishnan et al., 2021). Based on DCQ, in both mild and severe COVID-19 patients, a disease pathogenicity-specific pattern in neutrophils and CM was determined. The homeostatic association between the immune cells like T cells and DCs was dysregulated, further supported by the immunophenotyping of MNPs and its association with the T cells (CD4 + and CD8 + ). Due to the high molecular heterogeneity displayed by patients, we performed a sample stratification through SNF, based on both metabolomic and transcriptomic data. Using the personalized GSMM and SNF cluster-driven mild/moderate and severe COVID-19 groups, the FBA identified the transporter reactions that fuel TCA-cycle intermediates, such as the cis-aconitate, a-ketoglutarate, succinate, malate, and fumarate, that may play an essential role in COVID-19, which was further supported by sctMetabolomics of the monocytes. Gene and metabolite prioritization analysis indicated that monocarboxylate transporters (MCTs or SLC16 gene family) like SLC16A6 or nucleoside transporter genes like SLC29A1 that are expressed in the monocytes and metabolites such as a-ketoglutarate, succinate, malate, and butyrate could play a crucial role in COVID-19.
The DCQ predicted an increase in the neutrophils and CM in COVID-19 patients, which corroborated earlier findings (reviewed in Reusch et al., 2021). Following SARS-CoV-2 infection, an elevated number of both mature and immature neutrophils were observed in the nasopharyngeal epithelium (Chua et al., 2020), lung (Liao et al., 2020), and blood (Giamarellos-Bourboulis et al., 2020;Guan et al., 2020;Lourda et al., 2021). Elevation of neutrophil counts can lead to excessive formation of neutrophil extracellular traps (NETs), an extracellular network of chromatin, oxidant enzymes, and microbicidal proteins released by the neutrophils to control infection (Zuo et al., 2020), potentially contributing to hyperinflammation in COVID-19 patients. Moreover, a myeloid-driven atypical cytokine storm, mediated by the classical proinflammatory monocytes and neutrophils, could be the critical mediator of COVID-19 disease severity (Vanderbeke et al., 2021). The decrease in the number of circulating unconventional T cells, gd-T cells, and MAIT cells in severe COVID-19 patients compared with HC was in line with previous studies Parrot et al., 2020). Given that both cell types are important for antiviral responses, decreases in the blood could indicate their recruitment to tissue and increased cell death (Shi et al., 2021). Interestingly, the co-expression network of the cell-specific markers indicated a gradual loss (convalescent > mild > severe) of the association between marker genes of the T cells and DCs, as well as between marker genes of granulocytes and other immune cells at the system level. These findings are indicative of dysregulated crosstalk between T cells and DCs (Hivroz et al., 2012) and granulocytes with other immune cells , which could explain the insufficient control of the viral replication (Saichi et al., 2021) leading to severe COVID-19. Interestingly, the COVID-19 convalescent samples also showed the dysregulated crosstalk, further supporting a recent study where it has been shown that mild SARS-CoV-2 infection can leave an inflammatory imprint in the myeloid cells and can affect the macrophage effector functions (Bohnacker et al., 2022).
Several studies, including ours, have reported alterations in myeloid cells in blood and tissues of COVID-19 patients-with increased levels of CM, M-MDSC, and LDG (Cabrera et al., 2021) and a decrease in NCM, pDC, and DC1 cells-with these changes linked to platelet activation and degranulation, activation of the coagulation cascade, and overall disease severity (Agrati et al., 2020;Chilunda et al., 2021;Falck-Jones et al., 2021;Krishnan et al., 2021;Lourda et al., 2021;Schulte-Schrepping et al., 2020;Trombetta et al., 2021). In our study, using DCQ and immune profiling, we report similar changes in the MNP profile, and in addition to that, we also observed altered chemokine receptor expression in different MNP cells. Chemokines play a crucial role in modulating the immune landscape during viral infections by regulating innate and adaptive immune responses (Khalil et al., 2021). Systemic changes in chemokines' levels and receptors play critical roles in disease pathogenesis and progression (Majumdar and Murphy, 2020). Chemokine receptor signaling controls emigration from the bone marrow, cell survival, and migration into inflamed sites. Changes indicate the dynamic nature of an ongoing immune response to viral infection. Changes in MDSC cell subsets could have significant impacts on immune responses to SARS-CoV-2 (i.e., decreased T cell responses, antigen-presenting cell responses, etc.) (Emsen et al., 2022;Koushki et al., 2021). The fate of circulating monocytes during infection and inflammation is determined by differential expression of specific cell surface receptors, including CCR2, CX3CR1, and CCR5, which are known to modulate monocyte differentiation, migration, function, and survival (Chen et al., 2013b;Gschwandtner et al., 2019;K€ aufer et al., 2018;Mehta and Reilly, 2012). SARS-CoV-2 has been shown to upregulate CCR2, CCR5, and CCR1 expression in the human thoracic dorsal root ganglion (Khalil et al., 2021) and could directly or indirectly affect the levels of these receptors also on myeloid cells. In this study, we observed significantly increased levels of CCR2 and CCR5 on all monocyte subsets and in some DC subsets in the blood of COVID-19 patients, which are in line with the high levels of their ligands in the blood of COVID-19 patients Patterson et al., 2021a) and attraction to the inflamed lung. SARS-CoV-2-infected lung macrophages have also been shown to express the chemokine ligand 5 (CCL5), and the CCL5/CCR5 axis disruption by CCR5 blockade reduced inflammation, restored T cell lymphocytopenia, and reduced SARS-CoV-2 viremia in COVID-19 patients (Patterson et al., 2021b). Of note, the significant decrease in pDC frequency in the blood of COVID-19 patients and the downregulation of CX3CR1 in the remaining pDC in blood could possibly be explained by the binding of CX3CL1 (fractalkine) to CX3CR1. Fractalkine is a potent lung chemoattractant of immune cells and was previously found to be elevated in the serum of COVID-19 patients (Tong et al., 2020). Since DCs were reported in COVID-19 lungs (Liao et al., 2020), it is possible that DC disappearance from the circulation is related to their influx in lymph nodes and lung tissue to initiate and assist in the adaptive immune responses against SARS-CoV-2.
Another feature observed in COVID-19 disease is the altered plasma metabolic profile indicating a shift towards amino acid and fatty acid metabolism, energy, and lipid metabolism (Shen et al., 2020;Thomas et al., 2020;Wu et al., 2020). Previously, we and others reported that SARS-CoV-2 modulated the AKT/ mTOR/HIF-1 pathway that is a central regulator of glycolysis and other central metabolic pathways Chen et al., 2020;Codo et al., 2020). We have also shown that the inhibition of the glycolysis, glutaminolysis, and AKT/mTOR/ HIF-1 pathways inhibited the SARS-CoV-2 replication, as SARS-CoV-2 utilizes these pathways for its efficient replication Krishnan et al., 2021). The solute carrier (SLC) transporters are membrane-associated proteins that mediate the influx and efflux of metabolites (Schumann et al., 2020). Thus, SLCs play a significant role in maintaining glycolytic balance within a cell and significantly influence the development of metabolic diseases (Zhang et al., 2019b). In our previous study, immunophenotyping identified increased expression of the carbohydrate transporter GLUT1 (SLC2A1) in CD8 + T cells, IM, and NCM and amino acid transporter xCT (SLC7A11) in all monocytes in COVID-19 patients compared with HC, highlighting the role of metabolite transporters in COVID-19 disease (Krishnan et al., 2021). Moreover, the alterations of the transporters were more prominent in the monocytic population, which could potentially contribute to mediating the disease severity by modulating its impaired immune function through metabolic al-terations. In the present study, instead of explicitly looking into specific transporters, we tried to comprehensively identify the transporters uniquely associated with COVID-19. For this purpose, we have used metabolic flux-prediction models by integrating the blood transcriptomics and plasma metabolomics data and determined that SLC transporters, such as SLC29A1 (ENT1; equilibrative nucleoside transporter 1) and SLC16A6 (MCT7; monocarboxylate transporter 7), were uniquely essential in COVID-19 samples. The SLC29A1 is a nucleoside transporter and is critical for nucleotide biosynthesis. The SLC16 gene family, comprising 14 members of the MCT family, plays an important role in energy metabolism as it catalyzes the rapid transport of lactate and pyruvate across the cell membrane, which are essential components for glycolysis (reviewed in Halestrap, 2013). Although no study reported the role of the SLC transporter earlier, a co-expression analysis performed in microarray analysis data of influenza-infected pediatric patients reported upregulated expression of SLC16A6 (Zarei Ghobadi et al., 2019), suggesting its potential role in respiratory viral infections. The SLC16A6 has been reported to act as a transporter of ketone bodies like b-hydroxybutyrate out of the liver (Newman and Verdin, 2017) that is again internalized by tissues other than the liver and is utilized in the TCA cycle in mitochondria or fatty acid synthesis in the cytoplasm (Sheraj et al., 2021). A more recent study also reported SLC16A6 as a taurine transporter (Higuchi et al., 2022), indicating its role in energy metabolism. It is also possible that transporters such as SLC16A6, which was differentially upregulated in COVID-19, can also regulate the intracellular and extracellular levels of a-keto acids like a-ketoglutarate, which are essential components of the TCA cycle. We have also identified mitochondrial transporters such as SLC25A1 (CIC; citrate/isocitrate carrier) and SLC25A11 (OCG; mitochondrial 2-oxoglutarate/malate carrier) along with TCA-cycle intermediates such as a-ketoglutarate, citrate, succinate, fumarate, and malate differentially represented either in the cytosol or extracellular space. The prominent expression of SLC25A11 and SLC25A1 together with the upregulation of intracellular malate and a-ketoglutarate observed at the sctMetabolomic level is further supported by our GSMMs. Further, metabolite essentiality analysis identified butyrate, a microbiome-derived four-carbon short-chain fatty acid, as an essential metabolite in severe COVID-19, which plays an essential role in energy metabolism and intestinal homeostasis (Liu et al., 2018). Recently our integrative analysis of cytokines, metabolites, and microbiome features suggested a potential role of microbial-derived immunoregulatory processes in fatal outcomes in COVID-19 due to the failure of the negative feedback mechanism that should confine the cytokine storm (Albrich et al., 2022). Interestingly, MCTs are also involved in the butyrate transportation (Chang et al., 2014). The role of SLC transporters in the regulation of energy metabolism has been studied extensively for cancer but requires further attention in the context of viral infections that affect the metabolic pathways.
The system-wide metabolomics profile in COVID-19 patients with different disease states identified several biomarkers and shed mechanistic insight into SARS-CoV-2 infection. Even though metabolomics is closer to the phenotype, the metabolomics-only approach limits identifying the alterations of the metabolites or the metabolic pathways that can be influenced by gender, age, ethnicity, comorbidities, diet, etc. Application of GSMM by integrating the host transcriptomics and metabolomics and the virus biomass function can provide information about the altered metabolic reactions. Earlier seminal studies in the context of SARS-CoV-2, based on host-virus GSMM of human alveolar macrophages (Renz et al., 2020), human airway epithelial cells (Bannerman et al., 2021), and liver cancer cells (Huh7) (Yaneske et al., 2021), reported that the role of biosynthesis and degradation of amino acids, glycolytic shift, and energy metabolism are essential for progressive viral replication. As the SARS-CoV-2 was not detected in the blood cell populations in our cohort, we therefore hypothesized that the systemic metabolic alterations are most likely because of the bystander effect of the infection due to the inflammatory conditions. The specific metabolic environment may differ from the metabolic alteration during productive replication in the infecting cells.
In addition to using GSMMs for metabolic flux characterization and reporter metabolites of different COVID-19 subgroups, the generation of personalized metabolic models enabled the identification and prioritization of genes and metabolites at a patient-specific level in our study. We did a comprehensive characterization of COVID-19 at a personalized level integrating transcriptomic data combining information from biologically agreeable metabolic flux predictions and topological network analysis. Our data indicated a system-level altered flux in the mitochondrial TCA-cycle intermediates that can explain the dysregulated central metabolic pathways in COVID-19 and could potentially be linked with the glycolytic modulation by mTOR/ HIF-1 signaling and mitochondrial dysfunction (Krishnan et al., 2021). We also reported that by blocking two critical central metabolic pathways, glycolysis and glutaminolysis, viral replication and production were inhibited (Krishnan et al., 2021). Although systemic GSMM was not reported in other respiratory diseases caused by viruses, in our recent study to understand the natural control of HIV-1 infection, we observed the regulation of similar metabolic pathways, but the difference in the metabolic reactions potentiates a disease-specific contextualization of the metabolic flux (Ambikan et al., 2022). A viral disease-specific systemic GSMM atlas for other emerging and re-emerging viruses is underway.
The study has limitations that merit comments. Firstly, the RNA-seq data were generated from whole blood samples, which can cause bias in DGE analysis due to different blood cell populations. To overcome the bias, we computationally estimated the cell type proportions in each sample through DCQ. We used the information to adjust the bias to a minimum possible while performing DGE analysis. Second, although the DCQ analysis adapted from the original EPIC algorithm has worked here because of high sequence depth and sequence quality, it may not work in data with low sequence depth samples. It is recommended to cross-validate the finding with the original EPIC algorithm or by other phenotypic methods like flow cytometry. Third, the definition of hospitalized-mild and hospitalized-severe patient groups was based on oxygen consumption, and clinicians' definitions may leave a great deal of ambiguity as observed in the HCA. However, one of the strengths of our study was that we grouped the samples using high-throughput transcriptomics and metabolomics data by employing the multi-omics SNF methodology. The new patient groups were exclusively built on samples' multi-omics profiles, which brings high confidence than the clinical definition of the COVID-19 severity. Further, the group level GSMM and FBA were performed on the new patient groups (SNF) to minimize the ambiguity in the clinical definition, which may have caused erroneous results otherwise. Finally, the GSMMs was created by considering the whole blood as a single biological system. This poses an additional challenge as different cells perform metabolic functions in a slightly different manner. Ideally, the models should have been explicitly created for each cell type, but due to technological challenges, it was impossible to generate single-cell-type metabolomics data to accurately constrain the model's nutrient uptake of the system in various conditions, which is an essential factor in the procedure. Future developments enabling scRNA-seq integration in metabolic models to overcome the high sparsity that characterizes this data type will allow disentangling of the distinct metabolic flux distributions displayed by different cell types. However, to validate the findings observed, we developed and applied sctMetabolomics for the critical cell type, i.e., monocytes, further strengthening our GSMM observation.
In conclusion, with a combination of multi-modal systemswide transcriptomics, DCQ, and immunophenotyping of MNPs, we observed an essential role in coordinating the immune cells in COVID-19 severity. Novel data-driven patient stratification recapitulated many clinical properties previously known but further enabled us to uncover the potential mechanical consequences of COVID-19 infection in immune cells. The systems biology approaches, including network topology analysis and personalized GSMM, sctMetabolomics of monocytes, and gene mapping in published scRNA-seq data, displayed a progressive dysregulation of the central metabolic pathway (TCA cycle) concomitant with disease severity. Alterations in central carbon and energy metabolism, TCA-cycle intermediates like malate and a-ketoglutarate, and expression of metabolite transporters in monocytes were associated with disease severity at individual and group levels. These observations prompt the metabolic modulation of monocytes in specific COVID-19 patient groups and potentiate personalized targets for treatment in severe patients.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ujjwal Neogi (ujjwal.neogi@ki.se).

Materials availability
This study did not generate new unique reagents.

Data and code availability
All of the data generated or analyzed during this study are included in this published article and/or the supplementary materials. Created datasets and code are publicly available. The metabolomics data can be obtained from the https://dx.doi.org/10.6084/ m9.figshare.13336862. All the codes are available at https://dx.doi.org/10.5281/zenodo.6706499. Additional Supplemental Items are available from Mendeley Data at https://dx.doi.org/10.17632/5prg4tsz7n.1. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

Cohort description
In this study, we included two cohorts of hospitalized-mild (n=26) and severe (n=11) COVID-19 patients and COVID-19 PCR negative samples (n=31). The cohorts' clinical and demographic parameters were described elsewhere (Krishnan et al., 2021). The disease pathogenesis was defined based on oxygen consumption. The study was approved by the Regional Ethics Committee in Stockholm, Sweden, and performed in accordance with the Declaration of Helsinki.

METHOD DETAILS
Multi-omics analysis Whole blood was collected in Tempus TM blood RNA tubes. The RNA sequencing was performed using Illumina NovaSeq6000 in S4 mode. The transcriptomics data processing was performed as recently described by us . Differential gene expression analysis was performed using the R/Bioconductor package DESeq2 v1.26.0 (https://bioconductor.org/ packages/release/bioc/html/DESeq2.html). Confounding factors such as cell type proportions, BMI, and other possible hidden factors were adjusted while performing differential expression analysis. Bioconductor package RUVSeq v1.28.0 was used to compute the factor of unwanted variation and added to the DESeq2 design matrix (Risso et al., 2014). KEGG gene set enrichment analysis for differentially regulated genes was performed using PIANO v2.2.0 (V€ aremo et al., 2013) (nperm=500, geneset statistic=mean) and for communities using enrichr module of python package GSEAPY v0.9.16 (Chen et al., 2013a;Subramanian et al., 2005). Metabolomics and proteomics data generation were performed using the Metabolon HD4 (Metabolon Inc, NC, US) and Olink Immuno-Oncology Panel (Olink, Sweden), respectively (Krishnan et al., 2021). Normalized values of interleukins from Olink proteomics data are scaled down using the R function scale, and the mean values of each sample were visualized as boxplots. The single-cell RNAseq (scRNAseq) data was kindly obtained from Zhang et al. (2020), which had similar study groups to the present study.

Digital cell quantification (DCQ)
The deconvolution algorithm adapted from the method Estimating the Proportions of Immune and Cancer cells (EPIC) (Racle and Gfeller, 2020) was employed in the digital cell quantification (DCQ) procedure. The DCQ uses a reference cell type-specific gene expression profile and a list of signature genes specific for the cell types to calculate the cell type proportion. We have used RNA HPA blood cell gene data obtained from the Human Protein Atlas as a reference profile (Uhlen et al., 2017). The data consists of gene-level expression data of 18 blood cell types, and it is based on Human Protein Atlas version 20.1 and Ensembl version 92.38. Signature genes for the 18 blood cell types in the reference profile were downloaded from the online database CellMarker (Zhang et al., 2019a) and PanglaoDB (Franzé n et al., 2019). The transcript per million (TPM) transformed gene expression data of all genes from the samples were inputted to EPIC along with reference profile and signature gene list to estimate the cell proportion.