Human/SARS-CoV-2 genome-scale metabolic modeling to discover potential antiviral targets for COVID-19.

BACKGROUND: Coronavirus disease 2019 (COVID-19) has caused a substantial increase in mortality and economic and social disruption. The absence of US Food and Drug Administration-approved drugs for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) highlights the need for new therapeutic drugs to combat COVID-19.
METHODS: The present study proposed a fuzzy hierarchical optimization framework for identifying potential antiviral targets for COVID-19. The objectives in the decision-making problem were not only to evaluate the elimination of the virus growth, but also to minimize side effects causing treatment. The identified candidate targets could promote processes of drug discovery and development. SIGNIFICANT
FINDINGS: Our gene-centric method revealed that dihydroorotate dehydrogenase (DHODH) inhibition could reduce viral biomass growth and metabolic deviation by 99.4% and 65.6%, respectively, and increase cell viability by 70.4%. We also identified two-target combinations that could completely block viral biomass growth and more effectively prevent metabolic deviation. We also discovered that the inhibition of two antiviral metabolites, cytidine triphosphate (CTP) and uridine-5'-triphosphate (UTP), exhibits effects similar to those of molnupiravir, which is undergoing phase III clinical trials. Our predictions also indicate that CTP and UTP inhibition blocks viral RNA replication through a similar mechanism to that of molnupiravir.

Enzyme

Short/branched chain specific acyl-CoA dehydrogenase

Aldo-keto reductase family 1 member C4

Aquaporin-8

Aquaporin-9

CAD protein

Dihydroorotate dehydrogenase

Dihydropyrimidinase

Glycerol kinase

3-hydroxyisobutyrate dehydrogenase

Isocitrate dehydrogenase [NADP]

Glutamine-dependent NAD(+) synthetase

Ribose-phosphate pyrophosphokinase 3

Sodium/bile acid cotransporter

Excitatory amino acid transporter 4

Proton-coupled folate transporter

Sodium-dependent neutral amino acid transporter B(0)AT1

Sorbitol dehydrogenase

StAR-related lipid transfer protein 3

Thymidine kinase 2

Uridine-cytidine kinase 2

Uridine 5′-monophosphate synthase

Metabolite

3-Phospho-D-Glyceroyl Phosphate

2-Phospho-D-Glycerate

3-Hydroxyisovaleryl Coenzyme A

Acetyl Coenzyme A

2-Oxoglutarate

Alanine

Arginine

Arginyl-Alanyl-Alanine

Asparagine

Aspartic acid

Adenosine Triphosphate

N-Carbamoyl-L-Aspartate

Carbamoyl Phosphate

3,5-Diiodo-L-Thyronine 4-O-Sulfate

3(S)-Hydroxy-Tetracosa-6,9,12,15,18-All-Cis-Pentaenoyl Coenzyme A

10Z,13Z,16Z,19Z-Docosatetraenoyl Coenzyme A

Cytidine-5′-Triphosphate

Cysteine

Lauroyl Carnitine

(S)-Dihydroorotate

D-Glucose 6-Phosphate

D-Glucurono-6,3-Lactone

Glutamine

Glutamic acid

Glycine

Guanosine-5′-Triphosphate

Guanine

Homocarnosine

Propinol Adenylate

Histidine

Heptadecanoyl Coenzyme A

Isoleucine

Leucine

Leucyl-Aspartyl-Lysine

Lysine

Cholesterol-Ester-5,8,11-Eico

Methionine

Nicotinamide Ribotide

Orotate

Orotidine 5′-Phosphate

Phosphoenolpyruvate

Phenylalanine

Proline

5-Phospho-Alpha-D-Ribose 1-Diphosphate

Ubiquinone-10

Ubiquinol-10

Serine

Threonine

Tryptophan

Trans, Trans, Cis-Geranylgeranyl Diphosphate

Tyrosine

Uridine-5′-Monophosphate

Uridine-5′-TrIphosphate

Valine

Introduction

The coronavirus disease 2019 (COVID-19) pandemic was caused by the infectious spread of the novel new acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which belongs to as the genus β-coronavirus [1]. COVID-19 has resulted in a substantial increase in mortality and severe economic and social disruption worldwide [2]. According to the World Health Organization (WHO) COVID-19 dashboard [3], as of January 27, 2022, COVID-19 had been responsible for over 5.6 million deaths in 194 countries and over 360 million infections. Facing this unprecedented crisis, many research groups, industries, and governments have expended considerable effort and resources to develop vaccines and medications to combat COVID-19 [4]. Several prevention methods and treatments for COVID-19, such as mask wearing, social distancing, and vaccination, have been promulgated by the Taiwan Centers for Disease Control and Prevention and similar agencies in nearly every country [5]. Vaccination can prevent viral infection and reduce the spread of the disease. Authorized COVID-19 vaccines are now widely available. According to the WHO dashboard (https://covid19.who.int/), as of January 25, 2022, a total of over 9.6 billion vaccine doses had been administered. However, the absence of US Food and Drug Administration (FDA)-approved drugs against SARS-CoV-2 highlights an urgent need to design new drugs [6]. Several approaches [6], [7], [8], [9], [10], [11], [12] to drug screening and repurposing have been developed to identify potential agents for treating COVID-19.

The rapid identification of potential therapeutic targets for COVID-19 is essential. Computational methods and systems biology approaches can play key roles in the discovery of suitable drugs. Constraint-based modeling (CBM) has been successfully applied in fundamental research [13], [14], [15], the inference of oncogenes [16], [17], [18], [19], [20], [21], [22], the discovery of anticancer targets [23], [24], [25], [26] in oncology, microbial engineering [27], [28], [29], and other research fields. CBM uses data- and knowledge-driven constraints to identify feasible metabolic flux distributions for a given condition [13,14]. The SARS-CoV-2 Alpha variant has been incorporated into the human alveolar macrophage model iAB-AMØ-1410 [30], [31], [32], a genome-scale metabolic model (GSMM) of normal human bronchial epithelial cells [33], a human metabolic reaction model [34], and the human GSMM Recon 2.2 [35,36]. These metabolic models of infection have been applied to computationally identify targets for combating COVID-19. However, such numerical methods are inefficient for identifying target combinations. For example, Recon 2.2 accounts for 7785 reactions and 6047 species; thus, more than 30 million combinations must be analyzed to consider all two-target reactions through the numerical approach.

The present study integrated the SARS-CoV-2 Delta variant into the genome-scale human metabolic model Recon 3D [37] to present the viral infection. We developed a fuzzy multiobjective hieratical optimization framework based on a modification of Identifying Anticancer Target (IACT) framework [38,39] to mimic general wet-lab experiments and discover potential targets for treating COVID-19. The objectives were to limit viral biomass growth, maximize cell viability, and minimize metabolic deviation of the cells perturbed by potential treatments. The optimization method was employed in the infection model to consider target combinations of not only genes but also metabolites. Fuzzy set theory was applied to convert the multiobjective optimization problem into a maximizing decision-making problem, which was a mixed-integer trilevel optimization problem. Currently available commercial software tools cannot solve the fuzzy hierarchical optimization problem. Conventional genetic algorithms can be extended to solve the problem to obtain global solutions [40]. This study applied the nested hybrid differential evolution (NHDE) to solve the decision-making problem [38,39] to obtain optimal antiviral targets.

Materials and methods

Integration of SARS-CoV-2 into human metabolic model

Viruses are nonliving entities and do not have their own metabolism for reproduction; therefore, they are entirely dependent on the hosts that they infect. According to experimental studies, viral infection leads to significant metabolic alterations in the host, including increases in the glycolysis rate and changes in adenosine triphosphate (ATP) production [41], [42], [43], [44]. A key step in the viral replication is the synthesis of the viral biomass within the host cell, including structural proteins and genetic material.

Experimental observations have indicated that viral biomass synthesis causes significant metabolic flux changes in host cells and that metabolic perturbations can directly alter viral reproduction. Thus, GSMMs can be used to predict how various metabolic alterations affect viral reproduction and to identify antiviral targets to combat COVID-19. This study created an integrated host–virus (HV) GSMM to screen promising antiviral targets for COVID-19. The human GSMM Recon3D [37], which accounts for 5835 species, 10,600 reactions, and 2248 enzyme-encoding genes, was used as the host model (denoted as HT) in this study. Host cells infected with SARS-CoV-2 were integrated into the HV model. The gene and protein sequences of the SARS-CoV-2 B.1.617.2 Delta variant were downloaded from National Center for Biotechnology Information (NCBI) genome database [45] and used to generate the stoichiometric coefficients of a viral biomass objective function (VBOF). The VBOF is a pseudo-reaction simulating the production of virus particles, namely nucleotides and amino acids, and the associated energy metabolites required to produce virus particles. The stoichiometric information of viral lipids was not included in the VBOF because dynamic experimental data on viral envelopes are scarce.

The VBOF was generated in accordance with the seven steps described by Aller et al [44], Renz et al [30] and Delattre et al [35]. The pseudo-reaction occurs in cytoplasm and is expressed aswhere and are the stoichiometric coefficients of nucleotides (N) and amino acids (A), and calculated using protein and gene sequences, respectively. The stoichiometric coefficient of each metabolite in Eq. (1) is expressed as millimoles per gram of viral biomass and are calculated as follows: where the total molecules () of each nucleotide in a mole of virus particles is , where C is the genome copy number, is the nucleotide's frequency in the viral genome and is its frequency in the replication intermediates. The total molecules () of each amino acid per mole of virus particles is obtained using the protein sequence of structural () and non-structural () proteins and the formula , where and are the copy numbers of each structural and non-structural protein, respectively.

The amount of each nucleotide is converted into grams per mole of virus particles by using , where is the molecular weight (or molar mass, g/mol) of the nucleotide N. Similarly, the amount of each amino acid per mole of virus particles is converted into grams per mole of virus by , where is the molar mass of the amino acid A. The total mass of virus particles is calculated for the total of the genomic and proteomic mass by using .

Protein and RNA nucleotide polymerization require energy from ATP hydrolysis and pyrophosphate liberations. We revised the computation of the stoichiometric coefficient for the other components built from Aller et al [44], Renz et al [30] and Delattre et al [35] as follows. The total molecules of ATP () are calculated from the molecule of each structural () and non-structural () protein as follows. The polymerization of amino acid monomers requires approximately four ATP molecules per peptide bond, i.e. the constant k is defined as k = 4. Here, and are the length of amino acids for the k structural and non-structural protein, respectively. The copy numbers, and , are acquired from Delattre [35]. The stoichiometric coefficient of ATP is expressed millimoles per gram of virus: ATP is hydrolyzed into adenosine diphosphate (ADP), Pi, and H+ through the following reaction: ATP + H2O → ADP + Pi + H + free energy. The stoichiometric coefficients of ADP, Pi and H+ are equal to that of ATP: .

Amino acids are polymerized in cells to make polypeptides and proteins. The formation of peptide bonds consumes energy, which is obtained from the hydrolysis of ATP. Furthermore, amino acids polymerize through condensation polymerization; therefore, for every monomer added to a growing polymer chain, one molecule of water is also produced. As a result, the number of molecules of H2O should consider both the hydrolysis of ATP required for polymerization and water produced in the formation of the peptide bond. The stoichiometric coefficient of water can thus be calculated by the equationwhere represents the overall molecules of H2O required for ATP hydrolysis as from .

The polymerization of RNA nucleotide monomers to form the viral genome (+ssRNA or -ssRNA) releases PPi molecules (defined in the following expressions by the constant k =1). The number of molecules of PPi () required to form the viral genome (P) and replication intermediates (P) is calculated using the respective nucleotide counts: The frequencies, and of the virus genome and its replication intermediate are calculated from the viral RNA sequence. The stoichiometric coefficient of PPi in the VBOF is expressed in millimoles per gram of virus:

Antiviral target discovery problem

We developed a computer-aided strategy for screening potential therapeutic antiviral targets to combat COVID-19. The screening strategy was designed to identify not only antiviral enzymes, but also anti-metabolites. The antiviral target discovery (AVTD) framework described in Fig. 1 was formulated as a hierarchical optimization problem based on a modified form of IACT framework [38,39] to mimic a wet-lab experiment. The hierarchical screening procedure was formulated as a trilevel optimization problem consisting of an outer optimization problem with multiple objectives and subject to two loop inner optimization problems describing the characteristics of treated and perturbed cells.

Fig. 1

Flowchart of computer-aided screening for antiviral targets to combat COVID-19. (A) Gene and protein sequences of SARS-CoV-2 were downloaded from NCBI database. (B) A pseudo-reaction was constructed as a viral biomass objective function. (C) The human genome-scale metabolic network Recon3D was downloaded from Virtual Metabolic Human (https://www.vmh.life/). (D) The genome-scale metabolic model of the host-virus cells was created. (E) Flux distribution patterns for host cells were obtained from clinical data if available; otherwise, the template values were computed with flux balance analysis (FBA) and uniform flux distribution (UFD) problem without perturbation. (F) A set of antiviral targets was identified using the nest hybrid differential algorithm (NHDE), and used to compute the flux distributions of the host-virus cells during treatment. (G) The same targets were used to compute the flux distributions of perturbed host cells during treatment. (H) From fuzzy set theory, the flux distributions of the template, treated and perturbed cells are determined to evaluate the fitness of the targets. (I) If the fitness was unsatisfactory, the next antiviral targets were generated using the NHDE algorithm, and the procedure was repeated. (J) If the fitness was satisfactory, the target was identified as potential candidate.

The aim of the AVTD problem is to identify a set of potential therapeutic antiviral targets to remedy COVID-19. The identified antiviral targets require to fulfill three goals, i.e. the virus replication can be eliminated (Fig. 1F and H), the infected cells can restore to their healthy counterparts (Fig. 1E, F and H), and metabolic perturbation of the host cells are as small as possible to reach smaller side effects (Fig. 1E, G and H). The goals are formulated as fuzzy multiple objectives in the outer optimization problem as follows.

The first objective is the fuzzy minimization of the viral biomass, i.e. the viral biomass production () for treatment was achieved as close zero as possible. The second and third objectives are optimized the cell viability of the host-virus cells () and host cells () during treatment. The fourth and fifth objectives are used to optimize the flux distributions () and metabolite-flows () of the treated cells that achieved as similar the host template () as possible and the metabolic perturbations of host cells during treatment that got as close the template as possible, and are defined by fuzzy equality functions [46]. The inner optimization problems are expressed as follows.where v f/b is the forward-backward flux vector of reactions; N HV and N HT are the stoichiometric matrices for host-virus cells and host cells, respectively; and and are the positive lower bound (LB) and upper bound (UB) of the j forward-backward flux, respectively; and is the maximal cellular objective obtained through FBA. The approaches depend on the and of the i modulated reactions from the inner optimization problems. The AVTD framework can be employed to investigate reactions modulated on gene-centric or a metabolite-centric approach Fig. 2. illustrates the concepts underlying depicts the concept of both approaches, and a detailed example described in Additional File 2. In the gene-centric approach, the LBs and UBs of modulated reactions are restricted as follows:where Ω is the set of reactions regulated by isozymes determined using the gene-protein-reaction (GPR) associations, and modulation parameter δ is determined by a nested hybrid differential evolution (NHDE) algorithm. The metabolite-centric regulators modulate the synthesis reactions of the active metabolites. The LBs and UBs of modulated reactions for the i active metabolite are restricted as

Fig. 2

Illustrations of gene- and metabolite-centric approaches. (A) In the gene-centric approach, the reaction r1 is catalyzed by isozymes, E1 and E2, and the reaction r2 is regulated by E1. The isozymes, E1 and E2, are knocked out to inhibit r1, and r2 is also blocked by E1. (B) In the metabolite-centric approach, the metabolite, M1, is inhibited, and the synthesis reactions, r2, r3 and r4 are there by also inhibited.

Solving strategy

The outer optimization problem in the AVTD framework comprises various fuzzy objective functions. Following procedures similar to those used by Wang et al [38], the fuzzy multiobjective hierarchical optimization problem was converted into a maximizing decision-making problem using fuzzy set theory. The maximization problem is expressed as followswhere η and η are the membership grades of the VBOF, cell viability, and metabolic deviation, respectively, and w is a weighting factor. The transformation of each fuzzy objective function in the outer optimization problem describes as follows. A one-side linear membership function is applied to convert each fuzzy minimization objective in Eq. (13) into a decision criterion.where the LB and UB are provided by a user, and is the biomass growth rate of treated (TR) and perturbed (PB) cells calculated from the inner optimization problems. Fuzzy minimization of the VBOF is calculated using a method similar that in Eq. (18).

Converting each fuzzy maximization objective in Eq. (13) into a decision criterion requires a one-side linear membership function. The second and third goals in the outer optimization problem are combined to determine the cell viability of the treated and perturbed cells through the computation of mean-min operation [38] as Hence, overall cell viability is calculated as

A two-side linear membership function is employed to convert the fuzzy similarity objective from Eq. (13) into a decision criterion.

The decision criteria of fuzzy similarity for all fluxes in the GSMM are summed in all two-side membership functions as , where N is the total number of fluxes in the GSMM. The metabolite-flow alternations can be calculated using a similar equation: , where N is the total number of metabolites in the GSMM. The flow rate of the m metabolite is computed aswhere Ω is the set of species located in different compartments of the cell.

The metabolic deviations for treated and perturbed cells comprising flux and metabolite-flow alterations as

Hence, overall metabolic deviation is defined as

With the intersection of the membership functions, the fuzzy multi-objectives can transformed into the hierarchical fitness η by Eq. (17) for maximization. The optimality of the AVTD problem can be demonstrated using a procedure similar to that described by Wang et al [38]. The maximizing decision-making problem (17) is a mixed-integer trilevel optimization problem involving linear and quadratic programming. It is a high dimension NP-hard problem that no commercial software has solved. In this study, we modified the primal version of NHDE algorithm to solve the maximizing decision-making problem (17). The NHDE algorithm is a parallel direct search procedure and extended from the hybrid differential algorithm [47]. It has succeeded to solve the NP-hard problems such as strain design problems [29], oncogene inference problems [18], [19], [20], [21], [22] and anticancer target design problems [38,39]. The performance and solution quality of the NHDE algorithm depended on three key setting factors: tolerance ratio used in migration, population size, and maximum number of iterations. The tolerance ratio was set to 0.05. A population size of 50 was used and the default number of iterations was 100.

Results and discussions

Stoichiometric analysis

We used the stoichiometric coefficients of the amino acids and nucleotides in the viral pseudo-reaction for Alpha and Delta variants, and the reaction of cell growth of the host cell to compare the stoichiometry between infected and uninfected cells. Many of the amino acids and all of the nucleotides differed substantially in stoichiometry between the viral pseudo-reaction and biomass reaction of the host cells, as illustrated in Fig. 3 A. The coefficients of glutamic acid (Glu), histidine (His) and proline (Pro) for the infected cells decreased significantly. By contrast, the coefficients of tryptophan (Trp), asparagine (Asn) and tyrosine (Tyr) for the infected cells increased. The difference of stoichiometric coefficients for leucine (Leu), alanine (Ala) and glycine (Gly) were nonsignificant. The stoichiometric coefficients for the Delta variant and host cells are presented in Fig. 3B. The stoichiometric coefficients for Leu, Ala and Gly in both cells were higher than those of the amino acids and nucleotides (Fig. 3B). The stoichiometric coefficients of RNA nucleotides were more than 2.6-fold higher (log2 fold change ≥ 1.38). The stoichiometric coefficients for ATP and UTP were greater than those for guanosine-5′-triphosphate (GTP) and cytidine triphosphate (CTP), as illustrated in Fig. 3B.

Fig. 3

Stoichiometric analysis of amino acids and nucleotides in the viral biomass reaction. (A) Log2 fold changes of stoichiometric coefficients in amino acids and nucleotides of Delta and Alpha variants versus to host cells. (B) Stoichiometric coefficients of amino acids and nucleotides in the biomass reaction of Delta variant and host cells.

RNA and protein sequencing of the SARS-CoV-2 Alpha variant was also applied to determine the stoichiometric coefficients for the VBOF for comparison. The stoichiometry of the Alpha variant was nearly identical to that of the Delta variant, as illustrated in Fig. 3A. To compare the protein sequences of the viruses, we first used the Expasy Translate tool (https://www.expasy.org/resources/translate) to translate the RNA sequences into the corresponding protein sequences, which were identical to the protein sequences accessed from the NCBI database [45]. We then compared the variants’ sequence alignments, and identified few differences in structural and nonstructural proteins of these variations (Table 1 ). The results indicated that the structural spike protein (S protein) of both viruses consists of 1273 amino acids and exhibit nine variations, and the nucleocapsid has three variations. The nonstructural protein ORF1ab was found to comprise 7096 amino acids and exhibit seven variations. The altered locations of the proteins are listed in Table 1. We discovered the sequences of the structural membrane (M) and envelop (E) proteins and nonstructural proteins ORF7B and ORF10 to be completely identical.

Table 1

Protein sequences of SARS-CoV-2 Delta and Alpha variants. Amino acids in parentheses present the altered amino acids between the variants’ protein sequences. Gene and protein sequences were accessed from NCBI genome database (https://www.ncbi.nlm.nih.gov/nuccore/MZ724506).

Gene	Total No. of Amino Acids	Altered Location (Delta←→Alpha)
S	1273	142 (D←→G), 154 (K←→E), 218 (H←→Q),452 (R←→L), 484 (Q←→E), 614 (G←→D),681 (R←→P), 1071 (H←→Q), 1101 (D←→H)
N	419	3 (Y←→D), 204 (M←→R), 377 (Y←→D)
M	222	–
E	75	–
ORF1ab	7096	81 (Y←→H), 1567 (I←→T), 3646 (A←→T),4715 (L←→P), 5753 (I←→M), 6711 (R←→K), 6958 (R←→K)
ORF1a	4405	81 (Y←→H), 1567 (I←→T), 3646 (A←→T)
ORF8	121	69 (L←→S)
ORF7a	121	82 (A←→V)
ORF6	61	33 (T←→I)
ORF3a	275	26 (L←→S)
ORF7b	43	–
ORF10	38	–

Gene-centric approach

The NHDE algorithm [38,39] was applied to solve the maximizing decision problem in Eq. (17) by using the gene-centric approach defined in Eq. (15) to discover optimal enzyme targets. The algorithm was run several times to identify a set of single gene targets exhibiting the highest fitness among 2248 enzyme-encoding genes. The NHDE algorithm is a genetic algorithm that can obtain and rank targets by fitness. We identified five potential target enzymes with a membership grade for VBOF greater than 0.762, as listed in Table 2 . Two targets, dihydroorotate dehydrogenase (DHODH) and ribose-phosphate pyrophosphokinase 3 (PRPS1L1), participated in the metabolism on pyrimidine and purine, which are used for RNA and DNA synthesis. 3-hydroxyisobutyrate dehydrogenase (HIBADH) partially participates in the amino acid degradation related to protein synthesis. Aquaporin-9 (AQP9) and sodium/bile acid cotransporter (SLC10A1) participate in bile secretion. SARS-CoV-2 was recently discovered in the bile of a patient with severe COVID-19 [48,49]. The membership grade (η) for the VBOF reflects the percentage decrease in viral biomass production in treated cells. When η = 1, viral biomass production has completely stopped. On the basis of this computation, we determined that each of these single-target gene treatments could reduce viral biomass production by more than 76%.

Table 2

Optimal single antiviral target enzymes determined using the NHDE algorithm. η, and η are the membership grades for viral biomass growth rate, cell viability, and metabolic deviation, respectively, for the treated and perturbed cells. Higher η, and η indicate lower viral biomass growth rate, higher viability of treated and perturbed cells, and smaller metabolic alterations to the host, respectively. No. Drug denotes the number of drugs listed in the DrugBank database (https://go.drugbank.com/) that modulate the related gene.

Gene	Protein	ηVBOF	ηCV	ηMD	No. Drug	Participated Pathway
DHODH	Dihydroorotate dehydrogenase (quinone)	0.994	0.704	0.656	26	Pyrimidine metabolism
PRPS1L1	Ribose-phosphate pyrophosphokinase 3	1	1	0.247	0	Purine metabolism
HIBADH	3-hydroxyisobutyrate dehydrogenase	0.887	0.999	0.240	1	Valine, Leucine and Isoleucine Degradation
AQP9	Aquaporin-9	0.868	0.999	0.240	1	Bile Secretion
SLC10A1	Sodium/bile acid cotransporter	0.762	0.998	0.240	25	Synthesis of Bile Acids and Bile Salts

The enzyme DHODH catalyzes the oxidation of dihydroorotate (Dhor-S) to orotate (Orot) by using ubiquinone as an electron acceptor (Fig. 4 ). It is involved in producing pyrimidines, which are building blocks of DNA, RNA, and molecules such as ATP and GTP that serve as energy sources in cells. Downregulation of DHODH in the HV cells reduced the synthesis rate of Orot from 1.4 to 0.42 mmole/gDW/h. In addition, the flow rates of metabolites such as Dhor-S, Orot5p and Ump on the pyrimidine pathway decreased as illustrated in Fig. 4. The cell growth rate of the HV cells reduced by 99.4% after treatment with DHODH inhibitors. The cell viability grade η for treated and perturbed cells was 0.704, and the metabolic deviation grade η for the DHODH inhibition treatment was 0.656. The high metabolic deviation grades indicated that both flux patterns for treated HV cells and perturbed host cells were close to those of the host cells, thereby implying lower rates of side effect. Therefore, DHODH inhibitor may effectively combat COVID-19 while causing fewer side effect than other treatments cause (Table 2). Some studies have suggested that DHODH inhibitors can be used to treat autoimmune diseases [50] and cancers such as small-cell lung cancer [51] and acute myeloid leukemia [52]. Leflunomide is an approved DHODH inhibitor that widely used as a modest immune regulator to treat autoimmune diseases, in treating COVID-19 disease with a small-scale of patients [53]. An orally bioavailable compound, PTC299, is a potential anti-COVID-19 inhibitor of DHODH to reduce SARS-CoV-2 replication [54]. Based on the computational results, we observed DHODH to be a promising antiviral target for COVID-19; this result is consistent with previous reports [53], [54], [55], [56], [57]. The AVTD framework can be used to identify promising targets for treating COVID-19 and to decipher metabolic mechanisms of inhibition for treatment (Fig. 4). We searched DrugBank [58] and identified 26 drugs that inhibit DHODH, 25 that inhibit SLC10A1, and 1 that inhibits HIBADH and AQP9. These are all potential candidates for drug repurposing to combat COVID-19.

Fig. 4

Integration of a concise metabolic network regulated by the enzyme DHODH with viral replication. DHODH inhibition downregulates the conversion of Dhor-S to Orot in the host–virus cells. The numerical values in the box indicate the metabolite flow rates (mmol/gDW/h) for treated cells (TR), perturbed cells (PB), host–virus cells (HV), and host cells (HT).

Drug combinations can be used to increase therapeutic efficacy and reduce toxicity [59]. The use of such combinations may increase the success rate of drug repurposing. However, the wet-lab approach to identifying and validating effective combinations is limited by the excessive number of potential target combinations. Computer-aided screening may overcome the drawback at the expense of a large computational burden. This study employed the Recon3D GSMM [37] to integrate with the viral cell growth of SARS-Cov-2 and create an HV model accounting for 5835 species, 10,601 reactions, and 2248 enzyme-encoding genes. If the NHDE algorithm is applied for a single group of candidates to generate searching individuals, more than 2.5 million combinations must be evaluated. To identify two-target combinations and reduce the computational burden, the NHDE algorithm was applied for two groups of candidates to generate searching individuals. We performed a series of computations to obtain a set of two-target combinations and determine their optimal grades. We identified numerous combinations able to completely block viral biomass production (i.e., with a membership grades η = 1). However, the membership grades for cell viability and metabolic deviation varied. The 20 two-target combinations with the highest fitness values are presented in Fig. 5 . The membership grade for the VBOP for all treatments was equal to one (η = 1), and the combination of DHODH downregulation and thymidine kinase 2 (TK2) upregulation resulted in the membership grades for cell viability and metabolic deviation being 42% and 21% higher, respectively, than those achieved by single-target DHODH inhibition. TK2 is involved in the production and maintenance of mitochondrial DNA (mtDNA). Reduced mtDNA production can be caused by certain genetic variations. Upregulation of TK2 may prevent the dysregulation of biomass maintenance during treatment. As a result, targeting the combination of DHODH and TK2 achieved higher grades for cell viability and metabolic deviation. The 12 two-target combinations with the highest fitness in Fig. 5 resulted in improved metabolic deviation grades and thus may result in less severe side effects than other targets. The final eight combinations in Fig. 5 achieved maximal ATP production and minimal cell growth (η = 1) in the treated and perturbed cells but achieved metabolic deviation grades of approximately 0.35. Subsequently, we evaluated three-target combinations and determined that the three-target combination of DHODH, PRPS1L1, and PCYT2 increased the metabolic deviation grade to 0.771 but reduced the cell viability grade to 0.87.

Fig. 5

Membership grades for viral biomass objective function (VBOF), cell viability (CV), and metabolic deviation (MD) under various multi-target combination treatments.

Metabolite-centric approach

The viral pseudo-reaction in Eq. (1) used 20 amino acids and four nucleotides to produce viral biomass. Using the metabolite-centric approach, we first downregulated each one of these building metabolites to investigate the performance of each antimetabolite treatment, as illustrated in Fig. 6 . Downregulation of any of the amino acids or nucleotides could nearly terminate viral cell growth (η ≥ 0.91, except His) and achieve satisfactory membership grades for cell viability (η ≥ 0.67, except ATP). The membership grade for metabolic deviation reached 0.827 when aspartic acid (Asp) synthesis was downregulated (i.e., exhibited lower metabolic alterations compared to the template in host cells); as a result, Asp synthesis downregulation should result in fewer side effects than other treatments do. By contrast, downregulation of His was incapable of inhibiting viral biomass production, such that η = 0. Although the inhibition of ATP synthesis could reduce viral biomass production, it also prevented biomass maintenance in the normal cells, leading to membership grades for cell viability and metabolic deviation of 0.187 and 0.141, respectively.

Fig. 6

Membership grades for viral biomass objective function (VBOF), cell viability (CV), and metabolic deviation (MD) for treatments targeting each of 20 amino acids and four nucleotides.

Many antiviral drugs for treating SARS-CoV-2, such as molnupiravir and remdesivir, are undergoing phase III clinical trials. Both molnupiravir and remdesivir inhibit key enzymes of SARS-CoV-2, including viral RNA-dependent RNA polymerase (RdRp). Molnupiravir is a prodrug of synthetic nucleosides that mimic cytidine and uridine and exert antiviral action by introducing copying errors during viral RNA replication [60,61]. In this study, the inhibition of CTP and UTP exhibited effects similar to those of molnupiravir (Fig. 6). The prodrug remdesivir acts as an adenosine nucleoside triphosphate analog to interfere with the action of viral RdRp [62]. Our computational predictions indicate that ATP inhibition would exert an effect similar to that of remdesivir by blocking viral RNA replication but result in low cell viability and high metabolic alteration. We performed a series of computations to obtain a set of two-target combinations and their corresponding optimal grades (Fig. 7 ), and we compared the two-target combinations with their single-target counterparts (Fig. 6). His combined with Gly completely blocked viral biomass production (η = 1), and targeting the combination of ATP and GTP improved cell viability and attenuated metabolic deviation. Combinations involving one of six amino acids resulted in an enhanced metabolite flow rate and the combination of UTP and 2-phosphoglyceric acid (2pg) also improved cell viability and attenuated metabolic deviation.

Fig. 7

Membership grades for two-target combinations. The NHDE algorithm was applied to discover the most favorable two-target combinations, each of which consisted of one of 20 amino acids or one of four nucleotides with another metabolite in the model.

On the basis of the aforementioned computations, the building blocks of viral biomass—amino acids and nucleotides—were eliminated to diminish viral growth. In addition, we identified seven single antimetabolite targets that could nearly block viral biomass growth (Fig. 8 ) while achieving membership grades of cell viability and metabolic deviation greater than 0.67. These anti-metabolites related to metabolites in metabolic pathways regulated by DHODH are nucleotide derivatives that participate in the pyrimidine biosynthetic pathway to produce cytosine, thymine, and uracil for viral RNA replication, as illustrated in Fig. 4. Using these data, we identified various two-target antimetabolite combinations that obtained higher VBOF, cell vitality, and metabolic deviation grades than their single-target counterparts (Fig. 8).

Fig. 8

Membership grades for single-target antimetabolite and two-target combination treatments. These antimetabolites were determined by the NHDE algorithm and excluded the viral biomass building blocks—comprising 20 amino acids and four nucleotides—as candidates.

Conclusions

Authorized COVID-19 vaccines are now widely available. However, the absence of FDA-approved drugs against SARS-CoV-2 has highlighted an urgent need to manufacture and screen new drugs or repurpose existing drugs and identify promising targets for treating COVID-19. This study developed the AVTD optimization framework to mimic a wet-lab experiment for discovering antiviral gene and metabolite targets against COVID-19. The basic step in drug discovery and development processes is to discover potential antiviral targets. The identified antiviral enzymes and metabolites through AVTD framework can provide for biologists to carry on subsequent progressing in order to save a lot of times for screening procedures.

The gene and protein sequences of the SARS-CoV-2 Alpha and Delta variants were used to model the viral biomass growth reaction. Comparison of the protein sequences of the Alpha and Delta variants revealed few differences in their structural and nonstructural proteins, and the stoichiometries of the variants were nearly identical. We generated an integrated HV human GSMM using Recon3D and the viral stoichiometry. The AVTD framework employs a fuzzy hierarchical optimization method, which we used to identify potential antiviral targets using the integrated model. Potential treatments were evaluated in terms of their ability to limit viral biomass growth and metabolic deviation while maximizing cell viability in both treated and perturbed cells.

The AVTD framework identified not only gene regulator targets but also metabolite-centric targets. DHODH inhibitors reduced viral biomass growth by 99.4% and achieved cell viability and metabolic deviation grades of 70.4% and 65.6%, respectively. Consistent with our predictions, several articles have recently reported that DHODH inhibitors can potentially be used to treat COVID-19. We also identified two-target combinations that could completely block viral biomass growth while achieving even higher metabolic deviation grades. Some of the identified targets are modulated by drugs existing in the DrugBank database, as listed in Additional File 1. These drugs may serve as potential candidates for repurposing to develop new treatments for COVID-19.

Through a metabolite-centric approach, we identified 19 amino acids and four nucleotides that could nearly terminate viral biomass growth (except His) and achieve satisfactory membership grades for cell viability (except ATP). The inhibition of CTP and UTP exhibited effects similar to those of molnupiravir, which is undergoing phase III clinical trials for COVID-19 treatment. In addition, we determined that certain two-target combinations of antiviral enzymes and metabolites could achieve similar antiviral results with higher grades for cell viability and metabolic deviation.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.