Leveraging metabolic modeling to identify functional metabolic alterations associated with COVID-19 disease severity.

OBJECTIVE: Since the COVID-19 pandemic began in early 2020, SARS-CoV2 has claimed more than six million lives world-wide, with over 510 million cases to date. To reduce healthcare burden, we must investigate how to prevent non-acute disease from progressing to severe infection requiring hospitalization.
METHODS: To achieve this goal, we investigated metabolic signatures of both non-acute (out-patient) and severe (requiring hospitalization) COVID-19 samples by profiling the associated plasma metabolomes of 84 COVID-19 positive University of Virginia hospital patients. We utilized supervised and unsupervised machine learning and metabolic modeling approaches to identify key metabolic drivers that are predictive of COVID-19 disease severity. Using metabolic pathway enrichment analysis, we explored potential metabolic mechanisms that link these markers to disease progression.
RESULTS: Enriched metabolites associated with tryptophan in non-acute COVID-19 samples suggest mitigated innate immune system inflammatory response and immunopathology related lung damage prevention. Increased prevalence of histidine- and ketone-related metabolism in severe COVID-19 samples offers potential mechanistic insight to musculoskeletal degeneration-induced muscular weakness and host metabolism that has been hijacked by SARS-CoV2 infection to increase viral replication and invasion.
CONCLUSIONS: Our findings highlight the metabolic transition from an innate immune response coupled with inflammatory pathway inhibition in non-acute infection to rampant inflammation and associated metabolic systemic dysfunction in severe COVID-19.

Introduction

Since the COVID-19 pandemic began in early 2020, SARS-CoV2 has claimed more than six million lives world-wide, with over 510 million cases to date (World Health Organization, 2021). SARS-CoV2 infection can result in a range of symptoms from a cough, fatigue, nausea, and diarrhea, to more serious symptoms such as difficulty breathing, which can result in the need for oxygen supplementation or ventilation (CDC, 2022a). Long term side effects of SARS-CoV2 infection include chronic joint pain, continued difficulty breathing, depression, and in some cases permanent organ damage to the lungs, heart, and kidneys (CDC, 2022b). Because SARS-CoV2 is heavily glycosylated it is well equipped for immune evasion (Walls et al., 2016, 2020; Watanabe et al., 2020). Understanding how the host immune system responds differently between patients could provide an immune-directed treatment plan, as compared to many of the currently available treatments which focus on breaking the invasion-replication cycle.

Metabolism and the immune system are tightly linked, and by understanding the metabolic patient-to-patient differences in the context of SARS-CoV2 infection we can gain a better understanding of immune response (Alwarawrah et al., 2018). Therefore, to identify these potential therapeutic strategies, it is important to investigate how viral infection impacts host metabolism. Previous COVID-19 metabolomics research has focused on the impact of infection on host metabolism, as well as metabolic shifts associated with disease severity (Caterino et al., 2021; Páez-Franco et al., 2021; Sindelar et al., 2021; Stukalov et al., 2021; Thomas et al., 2020). As the pandemic continues, in order to reduce burden on the healthcare system we must not only investigate the physiological impact of infection on the body, but also elucidate the underlying functional metabolism that is causing these plasma-level differences in different patient responses in order to prevent hospitalization. Our work builds on previous investigations of metabolic shifts associated with COVID-19 disease severity (Hasan et al., 2021; Krishnan et al., 2021; Shen et al., 2020; Song et al., 2020; Su et al., 2020;Thomas et al., 2020).

To achieve this goal, we investigated metabolic signatures of both non-acute and severe COVID-19 samples by profiling the associated patient plasma metabolomes and building both severe and non-acute genome-scale metabolic models. By characterizing what functional metabolic pathways are upregulated in non-acute COVID-19 patients we hoped to identify protective metabolic pathways that assist in fighting infection, reducing immunopathology, and preventing disease progression. By identifying protective metabolic pathways and their associated pathways, this information can be leveraged as novel preventative therapeutic avenues. Furthermore, a better understanding of metabolic upregulation predictive of severe COVID-19 status can elucidate potential metabolic pathways whose inhibition could mitigate disease progression. Additionally, a deeper mechanistic understanding of SARS-CoV2 infection allows already available and approved treatment plans to be used more efficaciously.

We utilized both supervised and unsupervised machine learning approaches to identify key metabolic drivers that are predictive of COVID-19 disease severity and explored potential metabolic mechanisms that link these markers to disease progression. To more directly account for immune system response in disease progression, we have also analyzed the Interleukin 13 (IL-13) levels in the same patient cohort. IL-13 is a profibrotic cytokine that is involved in airway inflammation, and has recently been identified as a driver of COVID-19 severity (Donlan et al., 2021; Mohning et al., 2019). Our findings characterize the metabolic transition from an innate to adaptive immune response, the signatures of inhibited inflammatory pathways in non-acute COVID-19, and the metabolic byproducts of severe COVID-19 symptoms.

Materials and methods

Patients

Blood samples were collected in EDTA tubes from 84 adult patients who tested positive by PCR for SARS-CoV2 at the University of Virginia hospital between April and June 2020. Plasma prepared from the blood was stored at − 70 °C. A total of 48 of the samples were from out-patients and categorized as non-acute COVID-19, while 36 samples were categorized as severe COVID-19 based on the need for hospitalization and in some cases ICU and ventilator requirements (four and 25, respectively). All collected samples and subsequent data analysis had previous IRB approval.

Patient plasma preparation

Clinical information for the patients was obtained from the electronic medical records, and each patient was assigned a unique identifier for de-identification purposes. The protocol for collection of blood samples and de-identified patient information was approved by the University of Virginia Institutional Review Board (IRB-HSR #22,231 and #200,110). Metabolite profiling experiments were performed at Biomolecular Analysis Facility, University of Virginia School of Medicine. Briefly, plasma samples were thawed on ice and 50 µL of plasma was retained for the metabolome analysis. 200 µL of − 20 °C methanol was added to the plasma sample and shaken vigorously to inactivate any potential viruses. The samples were stored in − 80 °C immediately until extraction for metabolomics experiments. For extraction, 200 µL of − 20 °C methanol was added to each tube, vortexed and shaken vigorously for 30 min at 4 °C in a temperature controlled thermal shaker. Further 200 µL of chloroform and 400 µL of water were added, shaken vigorously and the top aqueous phase was recovered as a metabolite mixture of diverse chemical nature. Two aliquots of 350 µL each were created, one for the experiment and one as a backup for additional experiments. Each metabolite extract was dried overnight in speedVac and reconstituted in 60 µL of 0.1% formic acid in water. Prior to running the samples, the instrument was calibrated using Pierce FlexMix solution. Before running the actual experiment, LC–MS system was stabilized by running 3–4 wash runs followed by 4 blank runs to ensure stable background signal, and lastly a commercial amino acid mixture (Amino Acid Mixtures, 2022) was ran as a part of system stability analysis to ensure the column chromatograph was in good condition. Ten microliters from each tube were removed to create a pooled QC sample that was injected at the beginning and end of the MS sequence run and additional QC samples were injected after every 10 sample injections for continued assessment of chromatography quality.

Mass spectrometry preparation and analysis

For data acquisition, we employed a fully automated, AcquireX Intelligent Data Acquisition Workflow, specifically the deep scan workflow. First, an exclusion list is generated from a blank run. Then, an injection of the pooled QC sample followed by feature detection and component assembly populates the inclusion list with compounds detected in the samples. A series of iterative DDA injections follow. Each injection is informed from the previous one, minimizing redundant fragmentation spectra and maximizing relevant spectra and metabolite annotations.

Samples were analyzed in untargeted fashion by LC-HRMS. Samples were injected in randomized fashion via a Thermo Vanquish UHPLC and separation of the polar metabolites was achieved using Thermo Accucore C18 column (Thermo Scientific; 2.1 × 100 mm; 1.5 µm) maintained at 30 °C. The injection volume was 10 µL. For the 15-min gradient, the standard mobile phase for RPLC was A = 0.1% formic acid in water and B = 0.1% formic acid in methanol. The linear elution gradient was as follows: 0–8.0 min at 50% B, 8.0–13.0 min held at 98% B, 13.1 to 15.0 min revert to 0% B to re-equilibration for next injection at a flow rate of 0.25 mL/min. Spectra was acquired on Thermo IDX™ Tribrid MS, using both positive and negative mode. A heated electrospray ionization (HESI) source was operated at 3.5 kV and 2.5 kV for positive and negative modes, respectively. Ion source sheath gas was set at 35 and auxiliary gas at 7. Ion transfer tube temperature was maintained at 275 °C while vaporizer temperature was maintained at 320 °C. The instrument was set to acquire over the m/z range 67–1000, in full MS mode (1 µscan) at a resolution of 60,000 at a normalized AGC Target of 25% and 50 ms of maximum injection time was allowed. RF lens amplitude was set at 35%. Tandem MS/MS performed by applying quadrupole isolation with an isolation window of 1.6 m/z. Activation type was set at HCD and masses were fragmented with HCD Assisted Collision Energy (%) of 15,35,50. Fragment masses were detected by Orbitrap at a resolution of 15,000. A total of 140 × 2 sample runs (negative and positive mode), 30 QC runs (15 × 2 negative and positive mode), 8 ddMS2 runs (2 × 4) and 4 blanks (negative and positive) were collected.

LC–MS/MS data was processed using Compound Discoverer 3.1 (Thermo Fisher Scientific Inc, San Jose, USA) with the following settings. For precursor selection, MS(n-1) precursors were selected, and S/N Threshold was set at 1.5. Retention time alignment was performed using adaptive curve algorithm with maximum shift allowed for 2 min at a mass tolerance of 5 ppm. Compounds were detected at a mass tolerance of 5 ppm and minimum peak intensity threshold was set at 500,000. Preferred ions were set at either [M + H] + 1 or [M-H]-1, respectively, for positive or negative mode. Compounds were grouped at a mass tolerance of 5 ppm and RT tolerance of 0.2 min was allowed. A QC-based area correction was applied using a linear regression model with minimum allowable QC coverage of 50% and maximum QC area RSD allowed was 30%. Peak areas were normalized using constant median. Blank samples were used for detection and identification of background compounds and removed subsequently.

Compound annotations were performed by searching the ddMS2 masses in mzCloud. All metabolite identifications met MSI Level 2 requirements (Sumner et al., 2007). All spectra were matched with confidence using the HCD library from mzCloud using accurate matching of precursor ion and MS/MS spectra. Second search with either formula or exact mass was performed in ChemSpider, KEGG database and in-house database of IROA Mass Spectrometry Metabolite Library of Standards (MSMLS) (approximately, standards 550 compounds detected in both positive and negative mode under similar chromatographic setting). Statistical differential metabolite discovery was performed using healthy patients as a reference group and log2FC, p-value and adjusted p-value were computed after log10 transformation of the normalized peak areas using Compound Discoverer V3.1(Thermo Fisher Scientific, n.d.). If a metabolite was detected in both positive and negative modes, the positive mode was kept and the negative mode was discarded. Metabolites are considered changed if the p-value ≤ 0.05 and log2FC = 1 using. A total of 680 distinct metabolites was measured in each sample.

Biomarker identification

All metabolomics data were cleaned using median metabolite value filter followed by log transformation and autoscaling to normalize within and between samples. The Mann Whitney U test from the scipy python package was utilized to identify metabolites significantly enriched in either the non-acute or severe COVID-19 disease state (Nachar, 2008). Benjamini–Hochberg correction was used as the multiple test correction method from the statsmodel python package, in order to mitigate false discovery rate to less than 5% (Benjamini et al., 1995). From this metabolite list, 124 well-documented pharmaceutical or non-endogenous metabolites were removed to identify non-treatment based metabolic shifts. We used bootstrapping to account for the disbalance in groups. Using the supervised learning approach random forest, we assessed accuracy of the identified endogenous metabolites to discern between severe vs non-acute COVID-19 status using the holdout validation method (Breiman, 2001). Model accuracy was also assessed using receiver operating characteristic curve (ROC) and recall calculation.

Using the unsupervised learning approach Nonmetric Multidimensional Scaling (NMDS) and associated Euclidean distance, the previously identified endogenous metabolites were analyzed for significant COVID-19 status clustering using the vegan R package. A Permutational Multivariate Analysis of Variance (PERMANOVA) was utilized to identify if clustering of the severe vs. non-acute groups was significant again conducted using the vegan R package (Anderson, 2017).

In order to analyze IL-13 level correlation with COVID-19 disease severity, we subset the aforementioned endogenous differential metabolomics data to include the 25 patients with the highest IL-13 levels, as well as the 25 patients with the lowest IL-13 levels. An unpaired, one-tailed t-test was run in Excel on the select metabolites to assess significant differences based on IL-13 levels.

Pathway analysis

Pathway analysis was conducted to identify metabolic shifts significantly associated with disease state, as well as novel metabolic mechanisms underlying COVID-19 disease progression. The MetaboAnalyst version 5.0 toolbox via the pathway analysis module was utilized to run pathway enrichment analysis on the previously selected endogenous metabolites, enriched in either the severe or non-acute COVID-19 disease state (Xia et al., 2009). The platform utilizes the Kyoto Encyclopedia of Genes and Genomes (KEGG) to run over-representation analysis (ORA) to identify enriched metabolic pathways via a calculated enrichment ratio based on the actual versus predicted metabolite hits (Kanehisa and Goto, 2000). The significance of these calculated enrichment ratios is determined by an assigned p-value derived from a hypergeometric test using a binomial distribution.

Genome-scale metabolic modeling

To identify potential functional metabolic shifts that cannot be inferred from metabolomics alone, contextualized genome-scale metabolic models (GEMs) of a severe and non-acute COVID-19 disease state were constructed and analyzed using flux balance analysis (Edwards et al., 2002). The previously published whole-body GEM Recon3D was used as a base model for constructing a severe and non-acute COVID-19 disease model (Brunk et al., 2018). Recon3D is not cell-type specific, as the collected metabolomics data is from patient plasma samples and therefore represents non-cell-type specific metabolism. Recon3D's objective function is formulated to represent biomass maintenance, equivalent to homeostasis.

Prior to model pruning, the identified differential endogenous metabolites were manually matched to respective model metabolites. Metabolites not present within the model were not included. Exchange bounds of associated metabolites were adjusted to simulate open metabolic exchange (− 1000 lower bound, 1000 upper bound) for the respective disease state's differential metabolites.

Each newly contextualized Recon3D model was pruned using RIPTiDe. RIPTiDe prunes model reactions to only those required for a parsimonious flux solution. After pruning, FBA was run 500 times using Gapsplit to collect 500 flux samples from each pruned model because there is not a singular solution for satisfying the homeostasis objective function (Keaty et al., 2020). By taking 500 flux samples we are able to more accurately explore potential solution spaces. We compared conserved reaction flux values across the three models, using NMDS followed by a PERMANOVA, to determine if the three model's flux values were significantly separated. We then identified the top conserved reactions with predictive flux values capable of differentiating the severe versus non-acute models using random forest.

Results

An initial heatmap of the metabolomics data was constructed to compare severe vs. non-acute COVID-19 samples (Fig. 1). These results illustrate differential metabolic profiles associated with disease severity (Fig. 1). A total of 226 metabolites were identified as significantly differential after analyzing all metabolite values using the Mann–Whitney U-test to compare severe COVID-19 samples to non-acute COVID-19 samples. Once non-endogenous metabolites were removed, 80 metabolites were identified as significantly elevated in non-acute COVID-19 samples and 21 metabolites were identified as significantly elevated in severe COVID-19 samples. The remaining 101 metabolites were input into a random forest analysis and the resultant model was capable of predicting COVID-19 disease severity with 7.14% out-of-bag error rate. Compared to 14.19% out-of-bag error rate when all 680 metabolites were included, the removal of non-endogenous metabolites improved model classifier capabilities. NMDS and associated PERMANOVA showed significant grouping (R2 = 0.09, p-value < 0.001) of non-acute COVID-19 samples vs. severe COVID-19 samples, based on measured endogenous metabolite levels (Fig. 2) (Table 1).

Fig. 1

Heatmap of top 50 endogenous and non-endogenous differential metabolites for non-acute and severe COVID-19 patient plasma sample data

Fig. 2

a Top 20 endogenous differential features identified by random forest as most important for predicting COVID-19 disease severity. b Receiver operating curve for random forest generated COVID-19 status model predictability c) Non-metric Multi-dimensional Scaling (NMDS) based on all endogenous metabolites identified as significantly different between patient categories. (PERMANOVA: R2 = 0.09, p-value < 0.001)

Table 1

Sample metadata outlining patient medical background and COVID-19 related treatment

	Total	Average Age	Average BMI	Max Oxygen	Ventilator	Diabetes	Kidney Disease	Heart Disease	Lung Disease	Immunosuppression	Cancer	Remdesivir Treatment
Non-acute patient data
Female	27	56 (15–94)	33 (22–55)	10		11	5	7	1	2	2	5
Black	6	72	32	4		4	2	4
Asian	1	24		0
Other	10	47	32	10		3	3	1		2	1
White	10	58	35	8		4		2	1		1	5
Male	21	56 (33–84)	30 (14–62)	10		8	5	3	3	1	3	3
Black	6	64	30	10		3	1	2	2	1	1	1
Other	11	47	30	5		3	4	1			1	2
White	4	65	24	2		2			1		1
Severe patient data
Female	10	34 (21–55)	57 (28–75)	0	10	5	1	1		1	1	3
Black	1	55	60	0	1	1
Asian	1	33	65	0	1							1
Other	5	35	49	0	5	2					1	1
White	3	27	68	0	3	2	1	1		1		1
Male	26	30 (20–44)	61 (14–83)	15	20	12	4	4	3	2	1	4
Black	4	35	71	3	3	2	1		1		1
Other	10	28	46	15	8	5						2
White	12	30	70	4	9	5	3	4	2	2		2

Bold represents values for summary of sex, not subcalculated by race

To assess the interplay between IL-13 levels, COVID-19 disease severity, and functional metabolism, we did a high-level integration of IL-13 level data with the aforementioned metabolomics analysis. After sub-setting the metabolomics data based on the top 25 (average IL-13 pg/ml = 36.32) and bottom 25 (average IL-13 pg/ml = 5.26) IL-13 levels, unpaired one-tailed t-test analysis identified 14 significantly different metabolite levels from the 101 previously identified differential metabolites. From the 21 metabolites significantly elevated in severe COVID-19 samples, two metabolites were identified as significantly different when grouping by IL-13 levels, both of which are elevated in IL-13-high patients (Table 2). From the 80 metabolites identified as significantly elevated in non-acute COVID-19 samples, 11 metabolite levels were significantly altered when comparing high vs. low IL-13 levels (Table 2). All of these metabolites were significantly higher in low IL-13 level patients, except L-homocysteic acid, which was reduced in low IL-13 level patients.

Table 2

Metabolites significantly different (p-value < 0.05) when comparing high vs. low IL-13 patient metabolomics

COVID-19 Status	IL-13 Level	Enriched Metabolite
Severe	High	4-imidazolone-5-propanoate
		3-methylglutarylcarnitine
Non-Acute	High	L-homocysteic acid
	Low	(24R,24'R)-fucosterol epoxide alanyl-poly(glycerolphosphate)
		Erucamide
		L-Ascorbate 6-phosphate
		N2-acetyl-L-lysine
		11-Nitro-1-undecane
		β –leucine
		n-ribosylhistidine
		Trp-Phe
		trimethylsilyl N,O-bis(trimethylsilyl)serinate

After analyzing the differentially enriched metabolites associated with non-acute or severe COVID-19 using the KEGG database, several metabolic pathways appeared to be characteristic of the COVID-19 disease state. From the list of 80 metabolites that were differentially enriched in non-acute COVID-19 samples, higher levels of L-tryptophan, melatonin, 5-hydroxy-L-tryptophan, 3-hydroxyanthranilic acid, indoleacetaldehyde, and anthranilate are significantly associated with tryptophan metabolism (FDR = 0.07) (Fig. 3a).

Fig. 3

a Patient plasma levels of metabolites involved in tryptophan metabolism and melatonin synthesis (FDR = 0.07) b Patient plasma levels of the metabolite 5’-methylthioadenosine c Patient plasma levels of metabolites involved in ketone body biosynthesis d Patient plasma levels of metabolites involved in histidine degradation (FDR = 0.28) e Patient plasma levels of metabolites associated with shift in energy source

From the list of 20 metabolites that were differentially enriched in severe COVID-19, there is significant association between 4-imidazolone-5-propionic acid, imidazole-pyruvate, and methylimidazoleacetic acid with the histidine metabolism (FDR = 0.28) (Fig. 3b). Additional metabolites of interest were identified from the list of differential endogenous metabolites. Specifically, three acylcarnitines (hexanoylcarnitine, 3-methylglutarylcarnitine, and 9,12-hexadecadienoylcarnitine), L-gamma-glutamyl-L-leucine and D-galactonate in severe COVID-19 samples were identified as significantly higher compared to non-acute samples, whereas 5’-methylthioadenosine was identified as significantly higher in non-acute COVID-19 samples compared to severe cases (Fig. 4c).

Fig. 4

Metabolic pathway p-values are assigned based on how significantly the identified metabolites indicate pathway enrichment. a Metabolic pathways associated with non-severe COVID-19 metabolite predictors. Tryptophan metabolism and melatonin synthesis (FDR = 0.07). b Metabolic pathways associated with severe COVID-19 metabolite predictors. 19 disease sev (FD R = 0.28)

Genome-scale metabolic modeling provides scaffolding to investigate the intracellular metabolism that translates to extracellular metabolic byproducts which then result in physiological changes and disease progression. By curating a whole-body genome-scale metabolic model in the context of SARS-CoV-2 infection, using metabolomic data integration, we can better understand the functional metabolism and differential metabolic pathways that define severe and non-acute disease states. By accounting for metabolic inputs and byproducts using metabolomics data, we were able to construct a metabolic network that allows for insight to the functional and metabolic pathway shifts that differentiate disease state, resulting in changes to identified metabolite levels. Through this mechanistic understanding we gathered a snapshot of the disease's intracellular mechanisms, not just a picture of the physiological byproduct created from the metabolomics data. These two approaches working in tandem, metabolic modeling and metabolomics, create a more complete portrait of COVID-19's impact on the body from both an intracellular and extracellular perspective.

Of the original 21 metabolites identified has significantly higher in severe samples, ten of the metabolites were matched to annotated metabolites within the Recon3D model. Of the original 80 metabolites identified as significantly higher in the non-acute samples, 37 of the metabolites matched to annotated metabolites within the Recon3D model. After adjusting exchange bounds and pruning each model–severe COVID-19 and non-acute COVID-19–Recon3D's initial 10,650 reactions were reduced to 1119 reactions and 1066 reactions respectively for each contextualized model. The models share 878 reactions, with 241 reactions unique to the severe model, and 188 reactions unique to the non-acute model (Fig. 5a). NMDS and associated PERMANOVA showed significant model separation based on conserved reaction flux values (R2 = 0.22, p-value < 0.001) (Fig. 5b). Based on random forest, the top ten reactions for accurately classifying the severe versus non-acute models are: 5'-nucleotidase (IMP), sodium/ammonium proton antiporter, nucleoside-diphosphate kinase, 5'-nucleotidase (dAMP), transport of phenylalanine, asparagine, and methionine, adenosine facilitated transport in cytosol, fumarate:thiosulfate antiport (mitochondrial), and 5,10-Methylenetetrahydrofolate transport (cytosol to mitochondria), valine transaminase, transport of isoleucine, tyrosine, and tryptophan (Fig. 5c)..

Fig. 5

a Comparison of reaction presence between models, after pruning via Reaction Inclusion by Parsimony and Transcript Distribution b Non-metric Multi-dimensional Scaling comparison of median flux values across reactions conserved genome scale metabolic models (PERMANOVA: R.2 = 0.22, p-value < 0.001) c Heatmap of the median flux values, normalized within rows of the top reactions identified by random forest for differentiating severe vs nonacute models

These findings work concurrently with previous research showing positive correlation between IL-13 levels and fatty acid metabolism, as well as IL-13 level mediated COVID-19 disease severity (Knudsen et al., 2020; Zhu. et al., 2015). Within the analyzed patient cohort, increased IL-13 levels were positively associated with need for mechanical ventilation. Additionally, within COVID-19 mouse models, IL-13 neutralization increased survival rate and improved disease trajectory (Donlan et al., 2021).

Discussion

Metabolism plays an important role in COVID-19. While previous work has focused on using plasma metabolomics to identify the metabolic impact of COVID-19 on patient metabolism, as well as investigate disease progression, here we use plasma metabolomics data in conjunction with a genome-scale metabolic model to better analyze functional metabolic shifts associated with disease progression. Based on metabolomic analysis and computational metabolic modeling, we discovered that in non-acute COVID-19 patients' tryptophan, glutathione, pyrimidine and D-glutamine, D-glutamate metabolism are significantly enriched. These metabolic highlights may have physiological implications in regards to mitigated inflammatory response, decreased fibroblast accumulation, increased Vitamin-D levels, and cellular damage prevention. Conversely, in the case of severe COVID-19 patient data we discovered a significant enrichment in histidine (Fig. 4d). Physiologically these findings have been previously shown to correlate with increased viral titer and increased fatty acid oxidation(Doğan et al., 2021; Kimhofer et al., 2020; López-Hernández et al., 2021; Overmyer et al., 2021; Roberts et al., 2022; Thomas et al., 2020).

Non-acute COVID-19 metabolism

Differential enrichment of L-tryptophan, melatonin, 5-hydroxy-L-tryptophan, 3-hydroxyanthranilic acid, indoleacetaldehyde, and anthranilate in non-acute COVID-19 samples is significantly associated with tryptophan metabolism. Additionally, Trp-Phe is significantly enriched in both non-acute COVID-19 samples, as well as patients with low IL-13 levels. Tryptophan plays a role in neurotransmitter synthesis, anti-inflammatory pathways, as well as improved immune response via the gut microbiome (Gao et al., 2018; Krause et al., 2011; Wang et al., 2002). Increased tryptophan metabolism in non-acute COVID-19 patients presents a potential inhibitory mechanism for respiratory inflammation, a symptom associated with the cytokine cascade involved in SARS-Cov-2 infection (Dehhaghi et al., 2019). Notably, L-homocysteic acid is elevated in both non-acute COVID-19 samples, as well as IL-13-high samples. L-homocysteic acid is an endogenous neurotransmitter ligand of N-Methyl-D-aspartic acid, commonly known as NMDA. In regards to SARS-Cov2 infection, NMDA plays an integral role in neural plasticity, which could account for its high levels in the seemingly discordant categories of non-acute COVID-19 and high IL-13 levels (Jewett & Thapa, 2020). A potential mechanism is immune-mediated neural adaptation to disease pathology in order to avoid neurological symptoms that are sometimes reported in those infected with SARS-Cov2 (Marshall, 2020; Yirmiya & Goshen, 2011). In summary, increased tryptophan metabolism in non-acute COVID-19 could be directly upregulating anti-inflammatory pathways (Blasco et al., 2020).

Differential enrichment of 5′-methylthioadenosine (MTA) in non-acute COVID-19 samples could be a direct result of the body mitigating innate immune system inflammatory response. Mechanistically, MTA inhibits the pro-inflammatory cytokine TNF-α and increasing production of anti-inflammatory IL-10 (Hevia et al., 2004; Veal et al., 2004).

After investigating functional metabolic changes using a contextualized non-acute COVID-19 metabolic model three reactions stood out as carrying significantly higher flux compared to the severe COVID-19 metabolic model: nucleoside-diphosphate kinase, 5'-nucleotidase (IMP), and adenosine facilitated transport in cytosol. Nucleoside-diphosphate kinases (DNPK) plays an integral role in maintaining genomic stability, as they maintain the available pool of ribonucleotides and deoxyribonucleotides in the cell (Kapoor & Varshney, 2020). Previous research has shown that DNPK's can provide a protective effect in the case of cancer metastasis, a disease that also heavily relies on cellular replication (Lacombe et al., 2021). 5'-nucleotidase (IMP) dephosphorylates IMP and helps to maintain intracellular purine compound pools (Ipata & Tozzi, 2006; Pesi et al., 2021). Increased de novo purine synthesis has been linked to increased inflammatory response in COVID-19 patients. 5'nucleotidase catalyzes cytosolic purine degradation (Pieters & Veerman, 1988). Increased activity in 5'nucleotidase could be reducing cytosolic purine levels and thus reducing the cytokine response to infection (Xiao et al., 2021; Zhang et al., 2021). Lastly, adenosine is a signaling molecule which is present in higher concentrations during airway inflammation (Singh Patidar et al., 2018). Increased adenosine facilitated transport suggests that those who experience non-acute infection are removing excess adenosine and therefore reducing the subsequent inflammatory response to infection.

Severe COVID-19 metabolism

Significantly higher levels of imidazole-pyruvate, 4-Imidazolone-5-propanoate (also elevated in high IL-13 level samples), and methylimidazoleacetic acid in severe COVID-19 samples is associated with histidine metabolism, specifically histidine degradation (Fig. 4e). Additionally, 4-Imidazolone-5-propanoate is also elevated in high IL-13 level samples, indicating a potential immune system mediation of this metabolite. Histidine plays a key role in enzyme activation, specifically the activation of serine protease (Radisky et al., 2006). Host serine proteases have been shown to assist in viral entry into host cells, as well as viral spread of SARS-CoV2 (Seth et al., 2020). The serine protease alters the host S protein allowing for subsequent ACE2 receptor host-viral membrane fusion (Matsuyama et al., 2010). Increased histidine metabolism in severe COVID-19 samples could exemplify viral hijacking of host metabolism, as SARS-CoV2 continues to spread and uses host machinery against itself (Delattre et al., 2021).

Severe COVID-19 samples have significantly higher acetoacetate metabolite concentrations, which is associated with ketone body biosynthesis (Fig. 4d). When there is an increase in fatty acid oxidation there is a subsequent increase in ketone bodies produced by the liver (Botham Kathleen and Mayes, 2015). In addition, three acylcarnitines (hexanoylcarnitine, 3-methylglutarylcarnitine and 9,12-hexadecadienoylcarnitine) are significantly elevated in severe COVID-19 samples. 3-methylglutarylcarnitine is also elevated in high IL-13 level samples, connecting this metabolite to a potential inflammatory response pathway. Previous research has shown a positive correlation between increased IL-13 levels and increased fatty acid oxidation (NH et al., 2020). Acylcarnitines assist in the transport of fatty acids across the inner mitochondrial membrane for beta-oxidation. These two considerations in conjunction suggest an increase in fatty acid oxidation under severe COVID-19 conditions. Additionally, previous research has shown higher levels of fatty acids in SARS-CoV2 infected cells (Yan et al., 2019). However, it is important to note that increased levels of circulating acylcarnitines could be confounded by the significantly higher collective BMI of severe COVID-19 patients, as well as increased prevalence of heart disease (Kang et al., 2018; Mihalik et al., 2010; Ruiz et al., 2017; Smith et al., 2020).

Both L-gamma-glutamyl-L-leucine and D-galactonate are enriched in severe COVID-19 samples (Fig. 4f). L-gamma-glutamyl-L-leucine is associated with incomplete protein breakdown, and is excreted from the body via urine(Human Metabolome Database, n.d.; Jandke & Spiteller, 1986). D-galactonate is a metabolic byproduct of galactose breakdown, an energy-rich nutrient. An increase in galactose catabolism indicates a body that is in need of extra energy, while L-gamma-glutamyl-L-leucine’s association with incomplete protein breakdown could be a result of utilizing muscle for emergency energy. When considered in conjunction, these two metabolites suggest a metabolic energy shift in the severe disease state. Conversely, severe COVID-19 is associated with muscular atrophy, which could account for increased levels of L-gamma-glutamyl-L-leucine (Sagarra-Romero & Viñas-Barros, 2020). Additionally, at chronic high levels D-galactonate functions as a metabotoxin resulting in adverse health effects (Schlueter et al., 2018).

Our severe COVID-19 genome-scale metabolic model was able to capture functional metabolic shifts that were not identified from the metabolomics data alone. Specifically three reactions were highlighted via random forest: fumarate:thiosulfate antiport (mitochondrial), transport of phenylalanine, asparagine, and methionine, and transport of isoleucine, tryptophan, and tyrosine. High plasma levels of phenylalanine and tyrosine has been previously identified as biomarkers of increased disease severity in COVID-19 (Luporini et al., 2021; Shi et al., 2021). These findings indicate that our metabolic model is recapitulating previously identified markers of severe disease. Additionally, fumarate can function as a terminal electron acceptor in the electron transport chain when oxygen is not present (Spinelli et al., 2021). Increased fumarate:thiosulfate antiport in the mitochondria indicates reduced access to oxygen, which coincides with increased disease severity.

Conclusion

Enriched metabolites associated with tryptophan metabolism in non-acute COVID-19 samples, paint a metabolic picture of mitigated innate immune system inflammatory response, immunopathology-related lung damage prevention (Fig. 4a, b). Increased prevalence of histidine- and ketone-related metabolism in severe COVID-19 samples offers a potential mechanistic insight to musculoskeletal degeneration-induced muscular weakness, and host metabolism that has been hijacked by SARS-CoV2 infection to increase viral replication and invasion. Non-acute COVID-19 samples appear to have a reduced inflammatory response and a metabolic safety net to inhibit immunopathology-related lung damage, as opposed to severe COVID-19 samples whose metabolomes indicate take-over by the virus, resulting in a metabolic environment that is conducive to increased disease severity. These metabolic analyses and findings can be leveraged to improve and advance COVID-19 disease treatment, as well as lead to further investigation into metabolic predictors of disease advancement.