A computer-aided platform was developed in this study for screening potential therapeutic antiviral targets for treating heart cells infected with SARS-CoV2. This platform (Fig. 1) included a framework for cell-specific metabolic network reconstruction and AVTD. A reconstruction method, such as CORDA [36] or iMAT [37], can be used to reconstruct cell-specific GSMMs for HV and HT cells. To establish a generic VBR, the gene and protein sequence of SAR-CoV-2 from the NCBI database (https://www.ncbi.nlm.nih.gov/nuccore/) must be integrated with the generic human genome-scale metabolic network (GSMN) Recon3D to reconstruct cell-specific GSMMs. The RNA-seq expression of heart cells infected with SARS-CoV-2 is retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE150316) for reconstructing HV and HT models, respectively. Cell-specific GSMMs for HV and HT cells are then used for iteratively identifying antiviral targets on the AVTD platform.
2.1 Viral biomass reaction
The protein sequence of the SARS-CoV-2 Alpha (NC_045512) variant, which can be downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/nuccore/), can be used to generate the stoichiometric coefficients of protein and nucleotide polymerization in the VBR of this variant. Seven steps in the generation of protein and nucleotide polymerization of a VBR have discussed in [24, 30, 31]. Following such procedures, the gene and protein sequences of the SARS-CoV-2 are used to build the stoichiometric coefficients of RNA nucleotide and protein polymerization in the VBR. In this study, the ratio of the mass of lipids in the viral biomass to that in its host cell was used to estimate the stoichiometric coefficients of viral lipids in the VBR. The biomass reactions for HV and HT cells (cells infected with SARS-CoV-2 and healthy cells, respectively) are respectively expression as follows:
$$\begin{gathered} \sum\limits_{{i=1}}^{4} {S_{{{N_i}}}^{{HV}}{N_i}} +\sum\limits_{{j=1}}^{{20}} {S_{{{A_j}}}^{{HV}}{A_j}} +S_{{{H_2}O}}^{{HV}}{H_2}O+S_{{ATP}}^{{HV}}ATP+\alpha \sum\limits_{{p=1}}^{9} {S_{{{L_p}}}^{{HT}}} {L_p}\xrightarrow{{VBR}} \hfill \\ {\text{Virus-Biomass}}+S_{{ADP}}^{{HV}}ADP+S_{{Pi}}^{{HV}}Pi+S_{{PPi}}^{{HV}}PPi+S_{{{H^+}}}^{{HV}}{H^+} \hfill \\ \end{gathered}$$
1
where α is a mass ratio of lipids in the VBR relative to that of the host cells as follows:
$$\begin{gathered} \sum\limits_{{i=1}}^{4} {S_{{{N_i}}}^{{HT}}{N_i}} +\sum\limits_{{j=1}}^{{20}} {S_{{{A_j}}}^{{HT}}{A_j}} +S_{{{H_2}O}}^{{HT}}{H_2}O+S_{{ATP}}^{{HT}}ATP+\sum\limits_{p}^{P} {S_{{{L_p}}}^{{HT}}{L_p}} \xrightarrow{{VHOST}} \hfill \\ {\text{Host\_Biomass}}+S_{{ADP}}^{{HT}}ADP+S_{{{P_i}}}^{{HT}}{P_i}+S_{{P{P_i}}}^{{HT}}P{P_i}+S_{{{H^+}}}^{{HT}}{H^+} \hfill \\ \end{gathered}$$
2
where \(S_{{{N_i}}}^{{HV/HT}},S_{{{A_j}}}^{{HV/HT}},...,S_{{PPi}}^{{HV/HT}}{\text{ and }}S_{{{H^+}}}^{{HV/HT}}\)are the stoichiometric coefficients of nucleotides (Ni), amino acids (Aj), water (H2O), adenosine triphosphate (ATP), adenosine diphosphate (ADP), orthophosphate (Pi), diphosphate (PPi) and proton (H+) in the VBR of HV cells and the biomass reaction of a generic human GSMM of HT cells, respectively. The stoichiometric coefficients for each metabolite in the VBR can be calculated as follows:
$$\left\{ \begin{gathered} S_{{{N_i}}}^{{HV}}=1000\frac{{M_{{{N_i}}}^{{Tol}}}}{{M{W_{virus}}}} \hfill \\ S_{{{A_j}}}^{{HV}}=1000\frac{{M_{{{A_j}}}^{{Tol}}}}{{M{W_{virus}}}}\quad \hfill \\ S_{{{H_2}O}}^{{HV}}=1000\frac{{\left( {{k_{ATP}} - 1} \right)\left( {\sum {M_{{{A_j}}}^{{Tol}}} - 1} \right)}}{{M{W_{virus}}}} \hfill \\ S_{{ATP}}^{{HV}}=S_{{ADP}}^{{HV}}=S_{{Pi}}^{{HV}}=S_{{{H^+}}}^{{HV}}=1000\frac{{{k_{ATP}}\left( {\sum {M_{{{A_j}}}^{{Tol}}} - 1} \right)}}{{M{W_{virus}}}} \hfill \\ S_{{PPi}}^{{HV}}=1000\frac{{M_{{PPi}}^{{Tol}}}}{{M{W_{virus}}}} \hfill \\ \end{gathered} \right.$$
3
The total numbers of moles of the ith nucleotide, the jth amino acid, and PPi (\(M_{{{N_i}}}^{{Tol}}\), \(M_{{{A_j}}}^{{Tol}}\) and \(M_{{PPi}}^{{Tol}}\), respectively) are calculated from the corresponding molecule counts in the gene and protein sequences as follows:
$$\left\{ \begin{gathered} M_{{{N_i}}}^{{Tot}}={C_G}\left( {F_{{{N_i}}}^{G}+F_{{{N_i}}}^{R}} \right) \hfill \\ M_{{PPi}}^{{Tot}}={k_{PPi}}{C_G}\left( {\left( {\sum\limits_{i} {F_{{{N_i}}}^{G}} - 1} \right)+\left( {\sum\limits_{i} {F_{{{N_i}}}^{R}} - 1} \right)} \right) \hfill \\ M_{{{A_j}}}^{{Tot}}=\sum\limits_{k} {{C_{S{P_k}}}F_{{{A_j}}}^{{S{P_k}}}} +\sum\limits_{k} {{C_{N{P_k}}}F_{{{A_j}}}^{{N{P_k}}}} \hfill \\ \end{gathered} \right.$$
4
where the frequency \(F_{{{N_i}}}^{G}\) in the viral genome and the frequency \(F_{{{N_i}}}^{R}\) in the replication intermediate of each nucleotide Ni are calculated using the viral gene sequence retrieved from the NCBI database (https://www.ncbi.nlm.nih.gov/nuccore/); CG, \({C_{S{P_k}}}\) and \({C_{N{P_k}}}\) are the copy numbers of the gene sequence, the kth structural protein and the kth nonstructural protein, respectively. The frequency \(F_{{{A_j}}}^{{S{P_k}}}\) of amino acid Aj in the kth structural protein and the frequency \(F_{{{A_j}}}^{{N{P_k}}}\) of amino acid Aj in the kth non-structural protein are calculated using the viral protein sequence obtained from the NCBI database. The stoichiometric coefficient of water molecules is considered in protein polymerization [24, 29] is to account from the hydrolysis of ATP. However, water produces in the formation of the peptide bond during protein polymerization. Wang et al. [31] have revised the stoichiometric calculation to consider water produced in the formation of the peptide bond. Therefore, when Eq. (3) is used to obtain the stoichiometric coefficient of water, 1 M is deducted from the total number of moles of water produced during ATP hydrolysis.
The biomass reaction of HT cells [Eq. (2)] is used to calculate the mass ratio between the lipids and the biomass of HT cells as follows:
$$\begin{gathered} {\text{Ratio=}}\frac{{\sum\limits_{{p=1}}^{8} {S_{{{L_p}}}^{{HT}}{L_p}} }}{{{\text{Host_biomass}}}} \\ =\frac{{\sum\limits_{{p=1}}^{8} {S_{{{L_p}}}^{{HT}}{L_p}} }}{{\left( \begin{gathered} \sum\limits_{{i=1}}^{4} {S_{{{N_i}}}^{{HT}}{N_i}} +\sum\limits_{{j=1}}^{{20}} {S_{{{A_j}}}^{{HT}}{A_j}} +S_{{{H_2}O}}^{{HT}}{H_2}O+S_{{ATP}}^{{HT}}ATP{\text{+}}\sum\limits_{{p=1}}^{8} {S_{{{L_p}}}^{{HT}}{L_p}} - \hfill \\ \left( {S_{{ADP}}^{{HT}}ADP+S_{{{P_i}}}^{{HT}}{P_i}+S_{{P{P_i}}}^{{HT}}P{P_i}+S_{{{H^+}}}^{{HT}}{H^+}} \right) \hfill \\ \end{gathered} \right)}} \\ \end{gathered}$$
5
Here, Ni, Aj,…, and H+ denote as their corresponding molecular weights. Similarly, the mass ratio between the lipids and the viral biomass of HV cells is computed, and this ratio is assumed to be identical for HT and HV cells. Therefore, α is calculated as follows:
$$\alpha =\frac{{Ratio}}{{\left( {1 - Ratio} \right)}}\frac{{\left( \begin{gathered} \sum\limits_{{i=1}}^{4} {S_{{{N_i}}}^{{HV}}{N_i}} +\sum\limits_{{j=1}}^{{20}} {S_{{{A_j}}}^{{HV}}{A_j}} +S_{{{H_2}O}}^{{HV}}{H_2}O+S_{{ATP}}^{{HV}}ATP - \hfill \\ \left( {S_{{ADP}}^{{HV}}ADP+S_{{{P_i}}}^{{HV}}{P_i}+S_{{P{P_i}}}^{{HV}}PPi+S_{{{H^+}}}^{{HV}}{H^+}} \right) \hfill \\ \end{gathered} \right)}}{{\sum\limits_{{p=1}}^{8} {S_{{{L_p}}}^{{HT}}{L_p}} }}$$
6
In this study, the generic VBR in Eq. (1) was integrated with the generic human GSMN Recon3D to reconstruct cell-specific GSMMs for HV cells and HT cells, respectively. The reconstructed models were then used for iteratively identifying antiviral targets on the developed AVTD platform.
2.2 AVTD framework
Table 1
Hierarchical optimization framework based on four objectives for AVTD.
The objectives of the outer optimization problem are as follows: 1. To eliminate the viral biomass growth rate (VBGR) of HV cells as much as possible under the target treatment 2. To maximize the ATP production rate for treated HV cells and perturbed HT cells during treatment 3. To evaluate the fuzzy similarity between the metabolic patterns of treated HV cells and perturbed HT cells and those of HT cells 4. To evaluate the fuzzy dissimilarity between the metabolic patterns of treated HV cells and perturbed HT cells and those of HV cells The objectives of the inner optimization problem, which is subject to a constraint-based model, are as follows: 1. To conduct FBA and solve UFD problems for the treated HV cells 2. To conduct FBA and solve UFD problems for the perturbed HT cells |
We extended the AVTD framework developed in [31] to consider fuzzy dissimilarity objectives for evaluating the metabolic patterns of treated HV cells (denoted as TR cells) and perturbed HT cells (denoted as PH cells) with respect to those of HV cells. The AVTD framework described in Table 1 is aimed at mimicking a wet-lab experiment to identify targets for treatment. The four goals in the outer optimization problem are explained in the following part of this section. The first goal is to achieve the fuzzy minimization (\(\widetilde {{\hbox{min} }}\)) of the VBGR for TR cells, and this goal is expressed as follows:
$$\widetilde {{\mathop {\hbox{min} }\limits_{{\mathbf{z}}} }}\,v_{{VBGR}}^{{TR}} \approx 0$$
7
The second goal is to achieve the fuzzy maximization (\(\widetilde {{\hbox{max} }}\)) of the ATP production rate for TR and PH cells, and this goal is expressed as follows:
$$\left\{ \begin{gathered} \widetilde {{\mathop {\hbox{max} }\limits_{{\mathbf{z}}} }}\,v_{{ATP}}^{{TR}} \approx v_{{ATP}}^{{\hbox{max} }} \hfill \\ \widetilde {{\mathop {\hbox{max} }\limits_{{\mathbf{z}}} }}\,v_{{ATP}}^{{PH}} \approx v_{{ATP}}^{{\hbox{max} }} \hfill \\ \end{gathered} \right.$$
8
The third goal is to evaluate the fuzzy similarity (\(\widetilde {{{\text{similar}}}}\)) between the fluxes (vj) and metabolite flow rates (rm) of TR and PH cells and those of the HT template; this goal is expressed as follows:\(\left\{ \begin{gathered} \widetilde {{\mathop {{\text{similar}}\,}\limits_{{\mathbf{z}}} }}v_{j}^{{TR}} \approx v_{j}^{{HT}} \hfill \\ \widetilde {{\mathop {{\text{similar}}\,}\limits_{{\mathbf{z}}} }}r_{m}^{{TR}} \approx r_{m}^{{HT}} \hfill \\ \widetilde {{\mathop {{\text{similar}}\,}\limits_{{\mathbf{z}}} }}v_{j}^{{PH}} \approx v_{j}^{{HT}} \hfill \\ \widetilde {{\mathop {{\text{similar}}\,}\limits_{{\mathbf{z}}} }}r_{m}^{{PH}} \approx r_{m}^{{HT}} \hfill \\ \end{gathered} \right.\) (9)
The fourth goal is to evaluate fuzzy dissimilarity (\(\widetilde {{{\text{dissimilar}}}}\)) between the fluxes and metabolite flow rates of TR and PH cells and those of the HV template; this goal is expressed as follows:
$$\left\{ \begin{gathered} \widetilde {{\mathop {{\text{dissimilar}}\,}\limits_{{\mathbf{z}}} }}v_{j}^{{TR}} \approx v_{j}^{{HV}} \hfill \\ \widetilde {{\mathop {{\text{dissimilar}}\,}\limits_{{\mathbf{z}}} }}r_{m}^{{TR}} \approx r_{m}^{{HV}} \hfill \\ \widetilde {{\mathop {{\text{dissimilar}}\,}\limits_{{\mathbf{z}}} }}v_{j}^{{PH}} \approx v_{j}^{{HV}} \hfill \\ \widetilde {{\mathop {{\text{dissimilar}}\,}\limits_{{\mathbf{z}}} }}r_{m}^{{PH}} \approx r_{m}^{{HV}} \hfill \\ \end{gathered} \right.$$
10
In Eqs. (8)-(10), the decision variable z represents the gene-encoding enzyme determined by a nested hybrid differential evolution (NHDE) algorithm (Supplementary file 1). The fluxes (\(v_{j}^{{TR/PH}}\)) and metabolite flow rates (\(r_{m}^{{TR/PH}}\)) are to form the metabolic patterns of the TR and PH cells, and to obtain from the inner optimization problem by using each antiviral enzyme determined by the NHDE algorithm. The fluxes (\(v_{j}^{{HT/HV}}\)) and metabolite flow rates (\(r_{m}^{{HT/HV}}\)) of the HT and HV cell templates can be obtained from clinical experimental data (if available); however, genome-scale clinical data are currently unavailable. These templates can be obtained from HV and HT models, respectively, as standards for computing the inner optimization problems without regulation of an enzyme.
The flow rate of the mth metabolite of the TR and PH cells [Eqs. (9) and (10)] is computed as follows:
$${r_m}=\sum\limits_{{i \in {\Omega ^c}}} {\left( {\sum\limits_{{{N_{ij}}>0,j}} {{N_{ij}}{v_{f,j}} - } \sum\limits_{{{N_{ij}}<0,j}} {{N_{ij}}{v_{b,j}}} } \right)} ,m \in {\Omega ^m}$$
11
where Ωc represents the set of species located in various compartments of HT and HV cells and Nij is the stoichiometric coefficient of the ith metabolite in the jth reaction of a GSMM. The forward flux vf,j and backward flux vb,j of the jth reaction are calculated by applying FBA and UFD models in the inner optimization problem as follows:
$$\begin{gathered} \left\{ \begin{gathered} {\text{Treated HV model: }} \hfill \\ \left\{ \begin{gathered} {\text{FBA problem: }} \hfill \\ \mathop {\hbox{max} }\limits_{{{{\mathbf{v}}_{f/b}}}} {v_{BGR}} \hfill \\ {\text{subject to }} \hfill \\ {\text{ }}{{\mathbf{N}}^{HV}}\left( {{{\mathbf{v}}_f} - {{\mathbf{v}}_b}} \right)={\mathbf{0}} \hfill \\ \quad v_{{f/b,i}}^{{LB}} \leqslant {v_{f/b,i}} \leqslant v_{{f/b,i}}^{{UB}},\;i \notin {\Omega ^{TR}} \hfill \\ {\text{ }}v_{{f/b,j}}^{{LB,TR}} \leqslant {v_{f/b,j}} \leqslant v_{{f/b,j}}^{{UB,TR}},\;j \in {\Omega ^{TR}} \hfill \\ \end{gathered} \right\}\quad \left\{ \begin{gathered} {\text{UFD problem:}} \hfill \\ \mathop {\hbox{min} }\limits_{{{{\mathbf{v}}_{f/b}}}} \sum\limits_{{k \in {\Omega ^{Int}}}} {c_{k}^{{HV}}\left( {{v_{f,k}}+{v_{b,k}}} \right)} \hfill \\ {\text{subject to }} \hfill \\ {\text{ }}{{\mathbf{N}}^{HV}}\left( {{{\mathbf{v}}_f} - {{\mathbf{v}}_b}} \right)={\mathbf{0}} \hfill \\ {\text{ }}v_{{f/b,i}}^{{LB}} \leqslant {v_{f/b,i}} \leqslant v_{{f/b,i}}^{{UB}},\;i \notin {\Omega ^{TR}}{\text{ }} \hfill \\ \quad v_{{f/b,j}}^{{LB,TR}} \leqslant {v_{f/b,j}} \leqslant v_{{f/b,j}}^{{UB,TR}},\;j \in {\Omega ^{TR}} \hfill \\ {\text{ }}{v_{BGR}} \geqslant v_{{BGR}}^{*} \hfill \\ \end{gathered} \right\} \hfill \\ \end{gathered} \right\} \hfill \\ \left\{ \begin{gathered} {\text{Perturbed HT model:}} \hfill \\ \left\{ \begin{gathered} {\text{FBA problem: }} \hfill \\ \mathop {\hbox{max} }\limits_{{{{\mathbf{v}}_{f/b}}}} {v_{ATP}} \hfill \\ {\text{subject to }} \hfill \\ {\text{ }}{{\mathbf{N}}^{HT}}\left( {{{\mathbf{v}}_f} - {{\mathbf{v}}_b}} \right)={\mathbf{0}} \hfill \\ \quad v_{{f/b,i}}^{{LB}} \leqslant {v_{f/b,i}} \leqslant v_{{f/b,i}}^{{UB}},\;i \notin {\Omega ^{TR}} \hfill \\ {\text{ }}v_{{f/b,j}}^{{LB,TR}} \leqslant {v_{f/b,j}} \leqslant v_{{f/b,j}}^{{UB,TR}},\;j \in {\Omega ^{TR}} \hfill \\ \end{gathered} \right\}\quad \left\{ \begin{gathered} {\text{UFD problem:}} \hfill \\ \mathop {\hbox{min} }\limits_{{{{\mathbf{v}}_{f/b}}}} \sum\limits_{{k \in {\Omega ^{Int}}}} {c_{k}^{{HT}}\left( {{v_{f,k}}+{v_{b,k}}} \right)} \hfill \\ {\text{subject to }} \hfill \\ {\text{ }}{{\mathbf{N}}^{HT}}\left( {{{\mathbf{v}}_f} - {{\mathbf{v}}_b}} \right)={\mathbf{0}} \hfill \\ \quad v_{{f/b,i}}^{{LB}} \leqslant {v_{f/b,i}} \leqslant v_{{f/b,i}}^{{UB}},\;i \notin {\Omega ^{TR}} \hfill \\ {\text{ }}v_{{f/b,j}}^{{LB,TR}} \leqslant {v_{f/b,j}} \leqslant v_{{f/b,j}}^{{UB,TR}},\;j \in {\Omega ^{TR}} \hfill \\ {\text{ }}{v_{ATP}} \geqslant v_{{ATP}}^{*} \hfill \\ \end{gathered} \right\}{\text{ }} \hfill \\ \end{gathered} \right\} \hfill \\ \end{gathered}$$
12
where NHV and NHT are the stoichiometric matrices for the HV and HT model, respectively. These matrices and the corresponding gene-protein-reaction (GPR) association are reconstructed from Step A to Step D in Fig. 1. Moreover, \(v_{{f/b,j}}^{{LB}}\) and \(v_{{f/b,j}}^{{UB}}\) are the positive lower bound (LB) and positive upper bound (UB), respectively, of the jth forward flux and jth backward flux, respectively. The regulated LB and UB, namely, \(v_{{f/b,i}}^{{LB,TR}}\) and \(v_{{f/b,i}}^{{UB,TR}}\), respectively, depended on gene activation identified from the antiviral candidates generated by the NHDE algorithm [31, 38]. The regulation bounds based on GPR association can be expressed as follows:
$$\begin{gathered} {\text{Regulated bounds for }}{z_i}{\text{-th active gene in the GPR association:}} \hfill \\ {\text{Up-regulation:}} \hfill \\ \left\{ \begin{gathered} \left( {1 - \delta } \right)v_{{f,i}}^{{basal}}+\delta v_{{f,i}}^{{UB}} \leqslant {v_{f,i}} \leqslant \;v_{{f,i}}^{{UB}} \hfill \\ v_{{b,i}}^{{LB}} \leqslant {v_{b,i}} \leqslant \left( {1 - \delta } \right)v_{{b,i}}^{{basal}}+\delta v_{{b,i}}^{{LB}};\;{z_i} \in {\Omega ^{TR}} \hfill \\ \end{gathered} \right. \hfill \\ {\text{Down-regulation}}: \hfill \\ \left\{ \begin{gathered} v_{{f,i}}^{{LB}} \leqslant {v_{f,i}} \leqslant \left( {1 - \delta } \right)v_{{f,i}}^{{basal}}+\delta v_{{f,i}}^{{LB}} \hfill \\ \left( {1 - \delta } \right)v_{{b,i}}^{{basal}}+\delta v_{{b,i}}^{{UB}} \leqslant {v_{b,i}} \leqslant v_{{b,i}}^{{UB}};\;{z_i} \in {\Omega ^{TR}}{\text{\backslash }}{\Omega ^{IZ}} \hfill \\ \end{gathered} \right. \hfill \\ \left\{ \begin{gathered} \left( {1 - \varepsilon } \right)v_{{f,i}}^{{basal}} \leqslant {v_{f,i}} \leqslant \left( {1+\varepsilon } \right)v_{{f,i}}^{{basal}} \hfill \\ \left( {1 - \varepsilon } \right)v_{{b,i}}^{{basal}} \leqslant {v_{b,i}} \leqslant \left( {1+\varepsilon } \right)v_{{b,i}}^{{basal}};\;{z_i} \in {\Omega ^{TR}} \cap {\Omega ^{IZ}} \hfill \\ \end{gathered} \right. \hfill \\ {\text{Knockout}}:{\text{ }} \hfill \\ \left\{ \begin{gathered} {v_{f,i}}=0 \hfill \\ {v_{b,i}}=0;\;{z_i} \in {\Omega ^{TR}}{\text{\backslash }}{\Omega ^{IZ}} \hfill \\ \end{gathered} \right. \hfill \\ \left\{ \begin{gathered} \left( {1 - \varepsilon } \right)v_{{f,i}}^{{basal}} \leqslant {v_{f,i}} \leqslant \left( {1+\varepsilon } \right)v_{{f,i}}^{{basal}} \hfill \\ \left( {1 - \varepsilon } \right)v_{{b,i}}^{{basal}} \leqslant {v_{b,i}} \leqslant \left( {1+\varepsilon } \right)v_{{b,i}}^{{basal}};\;{z_i} \in {\Omega ^{TR}} \cap {\Omega ^{IZ}} \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered}$$
13
where \(v_{{f,i}}^{{basal}}\) and \(v_{{b,i}}^{{basal}}\) are the basal value of the ith forward and backward fluxes, respectively, obtained from the HV and HT templates; ΩIZ denotes the set of reactions regulated by isozymes identified using the GPR associations, and δ is the modulation parameter determined by the NHDE algorithm [31, 38]. The flux of a reaction catalyzed by isozymes remains around the basal level; thus, we set the flux ratio (ε) to 0.03 in the present study to restrict the flux value. The flux distributions and metabolite flow rates for the HV and HT cells can be used as HV and HT templates, respectively. These templates can be obtained by solving Eq. (12) without considering regulation of an enzyme.
In our previous study [31, 38], identical weighting factors (i.e., \(c_{k}^{{HV}}=c_{k}^{{HT}}=1\)) were used for solving UFD problems. In the present study, the RNA-seq expressions for HV and HT cells were used to not only reconstruct cell-specific GSMMs but also to set the weighting factors \(c_{k}^{{HV}}\) and \(c_{k}^{{HT}}\) for UFD problems to obtain uniform flux distributions. The weighting factors depended on quartile confidence classification using the RNA-seq expression of each cell. The four types of confidence reactions are as follows:
$$c_{k}^{{HV/HT}}=\left\{ \begin{gathered} \frac{1}{4},\;k \in {\text{high confidence}} \hfill \\ \frac{1}{2},\;k \in {\text{medium confidence}} \hfill \\ \frac{3}{4},\;k \in {\text{negativec confidence}} \hfill \\ {\text{1,}}\;k \in {\text{other confidence or non-gene-expression}} \hfill \\ \end{gathered} \right.$$
14
For a high confidence reaction, the smallest weighting factor is set to obtain a higher flux value in the UFD problem.
2.3 Maximizing Decision-Making Problem
The AVTD problem expressed in Eqs. (7)-(12) is a fuzzy multiobjective hierarchical optimization (FMHO) problem that can be transformed into a maximizing decision-making (MDM) problem by using fuzzy set theory to derive Pareto solutions (Fig. 2). One-side linear membership functions (blue dashed lines in Fig. 2) are defined to attribute fuzzy minimization and maximization and these functions are expressed as follows:
$$\begin{gathered} {\eta _{\hbox{min} }}=\left\{ \begin{gathered} 1,{\text{ if }}FV<LB \hfill \\ \frac{{UB - FV}}{{UB - LB}},{\text{ if }}LB \leqslant FV \leqslant UB \hfill \\ 0,{\text{ if }}FV>UB{\text{ }} \hfill \\ \end{gathered} \right. \hfill \\ {\eta _{\hbox{max} }}=\left\{ \begin{gathered} 0,{\text{ if }}FV<LB \hfill \\ \frac{{FV - LB}}{{UB - LB}},{\text{ if }}LB \leqslant FV \leqslant UB \hfill \\ 1,{\text{ if }}FV>UB{\text{ }} \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered}$$
15
where FV represents the fluxes or metabolite flow rates computed using the TR or PH model. The LB and UB are obtained using the corresponding HV or HT template; that is, LB = ST/4 and UB = 4ST. The term ST denotes the standard value for the HV or HT template used in the present study.
Two-sided linear membership functions are used to attribute fuzzy similarity (red line in Fig. 2) and fuzzy dissimilarity (green line in Fig. 2). The fuzzy similarity grade is derived using the equation as follows:
$$\begin{gathered} {\text{Left-hand side membership function:}} \hfill \\ \eta _{L}^{{TR/PH,HT}}=\left\{ \begin{gathered} 0,{\text{ if }}FV<LB \hfill \\ \frac{{FV - LB}}{{ST - LB}},{\text{ if }}LB \leqslant FV \leqslant ST \hfill \\ 1,{\text{ if }}FV=ST \hfill \\ \end{gathered} \right. \hfill \\ {\text{Right-hand side membership function:}} \hfill \\ \eta _{R}^{{TR/PH,HT}}=\left\{ \begin{gathered} 1,{\text{ if }}FV=ST \hfill \\ \frac{{UB - FV}}{{UB - ST}},{\text{ if }}ST \leqslant FV \leqslant UB \hfill \\ 0,{\text{ if }}FV>UB{\text{ }} \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered}$$
16
The fuzzy similarity grade \(\eta _{{MD}}^{{TR/PH,HT}}\) is obtained by combining the left-hand-side and right-hand-side membership functions as follows:
$$\eta _{{MD}}^{{TR/PH,HT}}=\hbox{max} \left\{ {\hbox{min} \left\{ {\eta _{L}^{{TR/PH,HT}},\eta _{R}^{{TR/PH,HT}},1} \right\},0} \right\}$$
17
.
The fuzzy dissimilarity grade is derived from the left-hand-side and right-hand side membership functions as follows:
$$\begin{gathered} {\text{Left-hand side membership function:}} \hfill \\ \eta _{L}^{{TR/PH,HV}}=\left\{ \begin{gathered} 1,{\text{ if }}FV<LB \hfill \\ \frac{{ST - FV}}{{ST - LB}},{\text{ if }}LB \leqslant FV \leqslant ST \hfill \\ 0,{\text{ if }}FV=ST \hfill \\ \end{gathered} \right. \hfill \\ {\text{Right-hand side membership function:}} \hfill \\ \eta _{R}^{{TR/PH,HV}}=\left\{ \begin{gathered} 0,{\text{ if }}FV=ST \hfill \\ \frac{{FV - ST}}{{UB - ST}},{\text{ if }}ST \leqslant FV \leqslant UB \hfill \\ 1,{\text{ if }}FV>UB{\text{ }} \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered}$$
18
The fuzzy dissimilarity grade \(\eta _{{MD}}^{{TR/PH,HV}}\) is obtained as follows:
$$\eta _{{MD}}^{{TR/PH,HV}}=\hbox{min} \left\{ {\hbox{max} \left\{ {\eta _{L}^{{TR/PH,HV}},\eta _{R}^{{TR/PH,HV}},0} \right\},1} \right\}$$
19
Eqs. (17) and (19) indicate that fuzzy dissimilarity is a complement of fuzzy similarity.
The AVTD problem is therefore transformed into an MDM problem by applying the membership functions as follows:
$$\left\{ \begin{gathered} \mathop {\hbox{max} }\limits_{{\mathbf{z}}} {\eta _D}=\mathop {\hbox{max} }\limits_{{\mathbf{z}}} {{\left( {\eta _{{CV}}^{{TR}}+\hbox{min} \left\{ {\eta _{{CV}}^{{TR}},\eta _{{CV}}^{{PH}},\eta _{{MD}}^{{TP}}} \right\}} \right)} \mathord{\left/ {\vphantom {{\left( {\eta _{{CV}}^{{TR}}+\hbox{min} \left\{ {\eta _{{CV}}^{{TR}},\eta _{{CV}}^{{PH}},\eta _{{MD}}^{{TP}}} \right\}} \right)} 2}} \right. \kern-0pt} 2} \hfill \\ {\text{subject to inner optimization problems}} \hfill \\ 1.{\text{ FBA and UFD problems for treated HV cells}} \hfill \\ {\text{2}}{\text{. FBA and UFD problems for perturbed HT cells}} \hfill \\ \end{gathered} \right.$$
20
where the decision objective ηD is a hierarchical criterion that the cell viability grade \(\eta _{{CV}}^{{TR}}\) of the TR model is used to achieve the first priority in the MDM problem. The cell viability grade \(\eta _{{CV}}^{{PH}}\) of the PH model and metabolic deviation grade \(\eta _{{MD}}^{{TP}}\) of the TR and PH models relative to their corresponding templates are considered the second priority of the decision objective. Membership grades for the MDM problem are defined as follows:
$$\eta _{{CV}}^{{TR}}={{\left( {\eta _{{VB{\text{R}}}}^{{TR}}+\hbox{min} \left\{ {\eta _{{VB{\text{R}}}}^{{TR}},\eta _{{ATP}}^{{TR}}} \right\}} \right)} \mathord{\left/ {\vphantom {{\left( {\eta _{{VB{\text{R}}}}^{{TR}}+\hbox{min} \left\{ {\eta _{{VB{\text{R}}}}^{{TR}},\eta _{{ATP}}^{{TR}}} \right\}} \right)} 2}} \right. \kern-0pt} 2}$$
21
$$\eta _{{CV}}^{{PH}}=\eta _{{ATP}}^{{PH}}$$
22
$$\eta _{{MD}}^{{TP}}=\frac{1}{2}\left( {\frac{{\left( {\eta _{{MD}}^{{TRHT}}+\eta _{{MD}}^{{PHHT}}+\eta _{{MD}}^{{TRHV}}+\eta _{{MD}}^{{PHHV}}} \right)}}{4}+\hbox{min} \left\{ {\eta _{{MD}}^{{TRHT}},\eta _{{MD}}^{{PHHT}},\eta _{{MD}}^{{TRHV}},\eta _{{MD}}^{{PHHV}}} \right\}} \right)$$
23
The membership grades \(\eta _{{VBR}}^{{TR}},\)\(\eta _{{ATP}}^{{TR}}\) and \(\eta _{{CV}}^{{PH}}\) in Eqs. (21) and (22) are obtained by using the one-side linear membership functions expressed in Eq. (15). The membership grades \(\eta _{{MD}}^{{TRHT}}\) and \(\eta _{{MD}}^{{PHHT}}\) are used to evaluate fuzzy similarity between the metabolic deviation for the flux patterns of TR and PH models and the HT template. The fluxes and metabolite flow rates of TR and PH models are used to compute the corresponding metabolic deviation grades according to the two-sided membership functions (Fig. 2) expressed in Eqs. (16) and (17). These grades are then used to compute overall metabolic deviation grades of the fuzzy similarity grades (\(\eta _{{MD}}^{{TRHT}}{\text{ and }}\eta _{{MD}}^{{PHHT}}\)). Similarly, the fuzzy dissimilarity grades (\(\eta _{{MD}}^{{TRHV}}{\text{ and }}\eta _{{MD}}^{{PHHV}}\)) between the metabolic deviation for the metabolic patterns of TR and PH models and the HV template are computed according to the two-sided membership functions (Fig. 2) in Eqs. (18) and (19), respectively.
The optimality and limitation of the transformation between the FMHO and MDM problems in Fig. 2 have been proved in a previous study [38]. According to the optimality theory, a Pareto solution of the FMHO problem can be obtained from the transformed MDM problem. The decision objective ηD [Eq. (20)] of the MDM problem is a hierarchical criterion. This criterion states that the cell viability grade \(\eta _{{CV}}^{{TR}}\) in Eq. (20) is the first priority in the MDM problem. Moreover, if the cell viability grade \(\eta _{{CV}}^{{PH}}\) or metabolic deviation grade \(\eta _{{MD}}^{{TP}}\) is less than \(\eta _{{CV}}^{{TR}}\), then one of the lowest grades in the set {\(\eta _{{CV}}^{{TR}}\), \(\eta _{{CV}}^{{PH}}\), \(\eta _{{MD}}^{{TP}}\)} becomes the second priority for decision-making. The MDM problem is a bilevel mixed-integer linear optimization problem and a high-dimensional NP-hard problem that cannot be solved using available commercial programs [39, 40]. We employed the NHDE algorithm to solve the aforementioned problem in this study. The NHDE algorithm is a parallel direct search algorithm that is an extended version of the hybrid differential evolution [41]. The implementation of the framework in this study are detailed in Supplementary file 1.