Protocol for condition-dependent metabolite yield prediction using the TRIMER pipeline

Summary This protocol explains the pipeline for condition-dependent metabolite yield prediction using Transcription Regulation Integrated with MEtabolic Regulation (TRIMER). TRIMER targets metabolic engineering applications via a hybrid model integrating transcription factor (TF)-gene regulatory network (TRN) with a Bayesian network (BN) inferred from transcriptomic expression data to effectively regulate metabolic reactions. For E. coli and yeast, TRIMER achieves reliable knockout phenotype and flux predictions from the deletion of one or more TFs at the genome scale. For complete details on the use and execution of this protocol, please refer to Niu et al. (2021).


SUMMARY
This protocol explains the pipeline for condition-dependent metabolite yield prediction using Transcription Regulation Integrated with MEtabolic Regulation (TRIMER). TRIMER targets metabolic engineering applications via a hybrid model integrating transcription factor (TF)-gene regulatory network (TRN) with a Bayesian network (BN) inferred from transcriptomic expression data to effectively regulate metabolic reactions. For E. coli and yeast, TRIMER achieves reliable knockout phenotype and flux predictions from the deletion of one or more TFs at the genome scale. For complete details on the use and execution of this protocol, please refer to Niu et al. (2021).

BEFORE YOU BEGIN
A metabolic network model is a mathematical abstraction of all known biochemical reactions and transmembrane transporters that occur within a living system for in silico computational predictions of metabolic dynamic behavior (Yurkovich and Palsson, 2015). Such networks provide a comprehensive view of all characterized metabolic processes by quantifying their metabolic reaction fluxes. Many methods for constraint-based network reconstruction and analysis can be used to simulate reaction fluxes through metabolic networks at the whole-genome scale (Bordbar et al., 2014). Among these, flux balance analysis (FBA) is a well-known technique that adopts linear programming to predict steady-state fluxes. For more accurate and robust prediction of target metabolic behavior under different conditions or contexts (e.g., for mutant strains due to gene deletions in metabolic engineering), these metabolic network models can also be integrated with a set of genetic regulatory rules, which can be modeled as a transcriptional regulatory network involving transcription factors (TFs) that may regulate metabolic reactions. Transcriptional regulation is often integrated via ''transcriptional regulatory constraints'' with various heuristics in metabolic network models for flux predictions (Covert and Palsson, 2003;Shlomi et al., 2007;Covert et al., 2008).
With the advent of high-throughput profiling technologies, genome-scale gene expression profiles can be easily obtained to help infer genetic regulatory rules. The recently developed hybrid model, Transcription Regulation Integrated with MEtabolic Regulation (TRIMER) (Niu et al., 2021), motivated by PROM (Chandrasekaran and Price, 2010), employs a Bayesian Network inferred from the gene expression data, built with the prior knowledge on regulatory relationships between TFs and genes to probabilistically model TF-gene regulatory rules. The resulting BN can be used to infer the joint probabilities in more general forms of Pr(gene(s)|TF(s)) than the marginal probabilities in forms of Pr(gene|TF) in PROM. These probabilities can subsequently be used to adjust the corresponding reaction flux constraints in the FBA formulation. This flexible probabilistic modeling of TF-gene regulatory relations enables modeling of the global effect of multiple TFs over their regulons, thereby better predicting the metabolic behaviors under different conditions, particularly when the focus is to predict fluxes with different TF Knockout (KO) strategies.
This protocol provides step-by-step explanations of applying TRIMER that integrates both TF-gene regulations and metabolic reaction regulations for in silico condition-dependent metabolic behavior modeling targeted for tasks such as knockout phenotype and knockout flux prediction. Specifically, we focus on analyzing genome-scale metabolic network models (e.g., IAF1260 for E. coli or IMM904 for S. cerevisiae) to predict metabolic fluxes and phenotypes of TF knockout strains for E. coli and yeast, to help identify knockout strains that overproduce the desired metabolites. The developed TRIMER pipeline can be extended for other organisms when both the corresponding metabolic reaction network model and transcriptomic expression data are available.
The main steps of the TRIMER pipeline comprise (1) learning a Bayesian network from gene expression data, (2) constructing flux constraints in predicting formulation, and (3) predicting metabolic fluxes. We will first present the necessary preparation steps and data requirements before applying TRIMER. Following this, we describe the critical steps for using TRIMER to predict metabolic fluxes under different TF KO conditions. Before introducing detailed implementation steps, we first list the corresponding input data, intermediate output files at each step, and final output results for TRIMER, as listed in Table 1, with corresponding examples in the TRIMER GitHub repository.

Installation of dependent software
Timing: 30-60 min

Preparation of metabolic model
Timing: 20-40 min Any organism with available gene expression data, a TF-gene interaction list and a metabolic model in the standardized data structure can be analyzed by TRIMER. For metabolic models, TRIMER requires them to be a MATLAB structure array and organized in the same way as required by the COBRA Toolbox, which is a popular module for constraint-based reconstruction and analysis of metabolic networks in MATLAB (Heirendt et al., 2019) (a typical metabolic model is shown in Figure 1). b. Field .grRules and .rules in the MATLAB array shown in Figure 1 correspond to the same set of GPR (Gene-Protein-Reaction) rules, written in two different forms. Illustrative examples are shown in Figure 2. Only GPR rules written in the way of the right column in Figure 2 can be processed by TRIMER. In field .rules, gene names are replaced by indices in the form of x(number), where the numbers represent the ordering of gene names contained in field .genes. In addition, logical operators are replaced by recognizable symbols in MATLAB. In case field .rules is missing, function cobra_to_trimer.m will try to create it based on field .genes and .grRules.
5. There are also many public databases of gene expression data and annotated TF-gene pairs.
When needed, interactions obtained from inference methods, such as GENIE3 (Huynh-Thu et al., 2010), TIGRESS (Haury et al., 2012), or Inferelator (Bonneau et al., 2006), can also be used to extend the interaction list. For this protocol, the authors use EcoMAC datasets (Carrera et al., 2014) at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4299492/. The TF-gene pairs in EcoMAC are obtained from the RegulonDBv8.1 database at http://regulondb.ccg.unam.mx/. Gene expression should be normalized and binarized, about which more details can be found in the step-by-step method details section. In TRIMER, exemplary binarized gene expression and TF-gene list can be found in the following two file paths: >TRIMER\source_data\Ecoli_model_EcoMac.mat The model is converted from the Boolean network model: iMC1010. The conversion code is in the path: The Boolean model iMC1010 is accessible at this link: https://systemsbiology.ucsd.edu/ InSilicoOrganisms/Ecoli/EcoliRegulations.

Configure mathematical programming solver
Timing: 1-3 min A mathematical programming solver should be specified and configured before metabolic analysis. Following TIGER (Jensen et al., 2011), TRIMER builds a customized MATLAB CMPI (Common Mathematical Programming Interface) for metabolic flux prediction based on the data structure detailed above. This CMPI defines a consistent and friendly interface for mathematical programming solvers, including GLPK and CPLEX for now.

MATERIALS AND EQUIPMENT
Hardware and operating systems: The authors used a Dell Inspiron 7370 computer with the Microsoft Windows 10 operating system. This computer has a 500 GB hard drive, an Intel Core i5-8250U CPU @ 1.60 GHz and 8 GB RAM. However, any reasonably up-to-date computer may be used to run all the code and any operating system can be used -Windows, Mac OS, or Unix/Linux. R software and required packages: The required packages should be installed for the first time when you run TRIMER. The authors used R (v4.0.3) and the following packages at the indicated versions when writing this protocol: bnlearn (v4.6.1) dplyr (v1.0.6) parallel (v4.0.3) purrr (v0.3.4) R.matlab (v3.6.2) binaryLogic (0.3.9) gRain (1.3.6) MATLAB software and required toolboxes: The authors used MATLAB 2016b and the following packages at the indicated versions when writing this protocol: COBRA (v2020) GLPKMEX (v4.62) Mathematical programming solver software: The authors used CPLEX version 12.8.0.

STEP-BY-STEP METHOD DETAILS
In this section, a comprehensive step-by-step protocol is laid out for applying TRIMER to a typical metabolic model of E.coli iAF1260 based on the EcoMac dataset. Interconnections of the main REAGENT Figure 3. All steps described are case-specific, but they can easily be adapted to any gene expression and TF-gene interaction dataset or metabolic model that the users wish to analyze. All the implemented MATLAB functions in TRIMER are well documented. Users can resort to MATLAB's help command for more details about function usages. Demo MATLAB code illustrating the complete procedure of applying TRIMER for biomass and indole flux prediction with TF-knockouts are accessible at the following paths: Construct TF-regulated gene network (TRN) Timing: 30-60 min We call the set of TF-gene regulations as TF-Regulated gene Network (TRN). In TRIMER, TRN is modeled by a Bayesian network, which is trained from pre-binarized data since TRN concerns 'on\off' states of TFs and genes.
Step 2   a. Quantile normalizes the raw expression data. b. Binarize the normalized data given a threshold. When writing this protocol, the authors chose the threshold to be 0.33 as suggested in the original PROM paper. In practice, the range of threshold is determined by evaluating PROM's performance. As PROM is computationally efficient and supported in TRIMER, the best threshold range for PROM can be chosen as the threshold for TRIMER. c. Perform Kolmogorov-Smirnov (KS) test over expression data for each TF-gene pair and prune out the ones that are tested to be independent. d. Here is an example of using data_preprocessing.m: CRITICAL: The binarization threshold is a very important hyperparameter, affecting TFgene pair selection and thereafter Bayesian network modeling of TRN.

(Under R environment) In TRIMER, bnlearn is used to learn Bayesian networks based on binary
expression data and TF-gene interaction lists. Network learning is comprised of two steps: structure learning and parameter fitting: a. Finding the global optimal structure by checking all possible directed acyclic graphs (DAGs) is computationally expensive as the cost grows exponentially with the number of nodes. Two structure learning strategies are implemented in TRIMER: Chow-Liu tree algorithm and Tabu search. i. Chow-Liu tree algorithm gives the optimal tree structure.
ii. Tabu search is a greedy searching method to learn general network structure. To further reduce the computation cost, TRIMER requires TF-gene pair as the searching space of Tabu search. b. Demo code for learning a Bayesian network using the EcoMac dataset is provided at the following path: c. The learned network is saved in the .bif format and will be used for further metabolic analysis in MATLAB. More details about .bif format can be found in the documents provided by the bnlearn package. A demo model is provided at the following path:

Infer conditional probabilities
Timing: 3-10 min Once we have learned a BN, we can make inference of all the relevant conditional probabilities which are in the form of Pr(gene(s)|TF(s)) and identified based on the TF-gene interaction list for one or Optional: TRIMER supports users to compute the conditional probabilities by relative frequencies of TF-gene state pairs as in PROM by calling the function prob_estimation.m.

Construct transcriptional constraints
Timing: 1-5 min With the list of affected genes and the associated regulatory values, new flux constraints can be derived for each TF KO condition. In this way, regulatory relationships between TFs and genes are integrated into predicting formulations (e.g., FBA).
4. Regulatory flux constraints are constructed by calling the function regulatory_bound.m. New flux constraints are obtained by multiplying the regulatory values computed from inferred conditional probabilities in the previous major step with the max/min flux values estimated via flux variability analysis (FVA) (Mahadevan and Schilling, 2003). For stability, flux constraints with values smaller than a threshold value (e.g., 10 -6 ) set by users are treated as zero. 5. Demo code for building new flux constraints is shown below: 6. TRIMER allows users to refine the regulatory flux bounds by calling the function regulatory_bound.m. Setting a minimum growth rate (biomass flux) requirement, this function returns the minimum set of regulatory flux constraints that cannot be satisfied when finding a feasible solution. This helps users to remove or adjust the bounds that over-constrained metabolic models. 7. Demo code for flux constraint refinement is shown below: >% Threshold value is set to be 1e-6 >% lb_est/ub_est -> constructed flux bounds >% vmax -> maximum/minimum flux values estimated by FVA > [lb_est,ub_est,vmax]=regulatory_bound(trimer,ko_tf,rxn_ affected_ko,rxn_prob_ko,'thresh',1e-6); 10. Besides the three formulations, TRIMER provides a set of useful functions that allow users to explore other formulations. Details about function usages can be found in their MATLAB help documents. Figure 6 shows the data structure used to represent all flux constraints. In this protocol, the authors take the implementation of sFBA as an example to show how to use these functions. The implementation can be achieved by the following steps: i. Replace the initial constraints with new regulatory constraints j. Make predictions for TF-knockout conditions: The data structures for constraint representations used in TRIMER are compatible with those used in TIGER. Flux constraints constructed by TRIMER can be directly added into the data structures in TIGER. This allows users to build the hybrid TRIMER-TIGER model where a part of flux constraints is constructed by TRIMER and the other part can be constructed by TIGER. This allows the modeling of the feedback effect from metabolites to genes. These feedback regulations cannot be obtained from gene expression data and are often expressed as Boolean rules. Demo code for using the TRIMER-TIGER model to predict flux phenotypes under various growth conditions can be found at the following path:

EXPECTED OUTCOMES
The main outcome of this protocol is to explain the procedure for using TRIMER to predict metabolic behaviors with TF knockouts. Here we present an example of the pipeline applied to the metabolic model for E. coli: iAF1260, following the workflow shown in Figure 3. The procedure begins with While this example is case-specific, a similar procedure described in this protocol can be applied to other organisms. Our code and step-by-step descriptions are intended to make the way of applying TRIMER more accessible to non-experts and to serve as a guide to other investigators for combining in silico flux simulations with regulatory network modeling.

QUANTIFICATION AND STATISTICAL ANALYSIS
To evaluate the performance of metabolic predictions by TRIMER, Pearson Correlation Coefficient (PCC) is used to evaluate the agreement between predicted fluxes and measured fluxes obtained experimentally. If TF knockout strains are tested under multiple growth conditions, we suggest using the normalized PCC, which is defined as the PCC between experimental flux ratios and predicted flux ratios. The flux ratio is computed by normalizing the knockout mutant strain fluxes by the corresponding wild-type flux in the corresponding growth condition. In this way, the normalized PCC removes the influence due to different growth conditions and can better illustrate how well is the knockout flux predictions ( Table 2 show an example of predicted fluxes and experimentally measured ground-truth fluxes with the corresponding PCCs and normalized PCCs). >TRIMER\grRate_KO.m.

LIMITATIONS
In TRIMER, the set of affected genes for a set of TFs considered is determined by the annotated TFgene pair list. In this protocol, RegulonDB, a gold-standard TF-gene list with strong experimental validation support is used. Such a validated list is important for achieving reasonable prediction performances. However, it should be noted that the gold-standard interactions are not necessarily perfect and may contain false-positive interactions. Similarly, the low-confidence interactions identified from other sources, for example by computational inference methods, could be either false positives or true interactions that have not been validated yet (Chandrasekaran and Price, 2013). This means that potential performance improvement can be still possible with better interaction lists (Wang et al., 2017). However, currently there is no commonly-adopted procedure to refine TF-gene interactions.
For metabolic network modeling, there is no standard operating procedure for determining the uptake rates of nutrients represented by lower bounds of the corresponding reactions. In the absence of in vivo uptake rates obtained from time-course metabolomic experiments, it is suggested to use uptake rates measured under wild-type conditions, which may not achieve the best flux predictions.

TROUBLESHOOTING Problem 1
Error occurs when using cobra_to_trimer.m to processGPR rules.

Potential solution
Only when GPR rules are in the standard form shown on the left column of Figure 2, running co-bra_to_trimer.m is guaranteed to transform them into the MATLAB recognizable form as shown on the right column of Figure 2. Otherwise, users are suggested to manually transform them into the standard form first before running cobra_to_trimer.m.

Problem 2
Failure to call the R script for probability inference in the MATLAB environment. The unit of fluxes is mmol/gDCW/hr; fluxes are measured under aerobic and anaerobic conditions with glucose as the only substrate (Covert et al., 2004); uptake rates of glucose and oxygen were set to be 8.5 and 14.6 mmol/gDCW/hr when making predictions for the aerobic condition; and 20.8 and 0 14.6 mmol/gDCW/hr for the anaerobic condition, respectively.

Potential solution
It is likely that the input R environment path is incorrect. Users should specify the installation path of the R environment when using TRIMER. The installation path of the R environment varies depending on the user's computer setup. Here is an example to define the path:

Problem 3
The running time and required memory to compute flux bound constraints for affected reactions is high.
Potential solution TRIMER supports two ways of computing the flux bounds. The first relies on estimating probabilities in the form of P(gene=1|TF=0) for each TF-gene pair, which is computationally efficient. In the second way, the probabilities are estimated in the form of P(genes=p|TFs=0) and are inferred by considering each valid gene-TF state pair denoted by p. It is possible that for a given reaction, the number of genes that affect this reaction and are also affected by the TF knockout is large. This means that there are an exponentially growing number of probability values to estimate with respect to the number of genes. Thus, implementing the second way can be computationally expensive. To avoid this computational problem in TRIMER, the default setup to infer conditional probabilities is to estimate P(gene=1|TF=0) for each TF-gene pair when the number of genes is higher than 12. Users can change this setting based on their need and computer setups.

Problem 4
There may be two potential issues when running functions for inferring conditional probabilities: (1) by multiple function calls for inference, the returned results are different; (2) the returned probability values may not satisfy the normalization axiom.

Potential solution
The default inference method in bnlearn is by logic sampling, an approximate inference algorithm, so that the returned results by multiple separate function calls can be different and they may not add up to 1. If these are observed, the user should increase the sample size of logical sampling to achieve more robust estimates. The default choice of the sample size is 5000 x log 10 (number of BN model parameters).
Problem 5 MATLAB crashes when using the ROOM formulation to predict fluxes for some knockout conditions.
Potential solution MATLAB often freezes when mathematical programming solvers (CPLEX or GLPK) are unable to reach optimality, and instead iterate the same solution repeatedly. It is a common phenomenon when solving mix-integer programming problems, such as ROOM. To address this issue, users can change the settings of the solvers to avoid this. For example, users can set the maximum running time of the solver to 10 min.

RESOURCE AVAILABILITY
Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Xiaoning Qian, xqian@ece.tamu.edu.

Materials availability
This study did not generate new unique reagents.