Genome-scale metabolic modeling of responses to polymyxins in Pseudomonas aeruginosa

Abstract Background Pseudomonas aeruginosa often causes multidrug-resistant infections in immunocompromised patients, and polymyxins are often used as the last-line therapy. Alarmingly, resistance to polymyxins has been increasingly reported worldwide recently. To rescue this last-resort class of antibiotics, it is necessary to systematically understand how P. aeruginosa alters its metabolism in response to polymyxin treatment, thereby facilitating the development of effective therapies. To this end, a genome-scale metabolic model (GSMM) was used to analyze bacterial metabolic changes at the systems level. Findings A high-quality GSMM iPAO1 was constructed for P. aeruginosa PAO1 for antimicrobial pharmacological research. Model iPAO1 encompasses an additional periplasmic compartment and contains 3022 metabolites, 4265 reactions, and 1458 genes in total. Growth prediction on 190 carbon and 95 nitrogen sources achieved an accuracy of 89.1%, outperforming all reported P. aeruginosa models. Notably, prediction of the essential genes for growth achieved a high accuracy of 87.9%. Metabolic simulation showed that lipid A modifications associated with polymyxin resistance exert a limited impact on bacterial growth and metabolism but remarkably change the physiochemical properties of the outer membrane. Modeling with transcriptomics constraints revealed a broad range of metabolic responses to polymyxin treatment, including reduced biomass synthesis, upregulated amino acid catabolism, induced flux through the tricarboxylic acid cycle, and increased redox turnover. Conclusions Overall, iPAO1 represents the most comprehensive GSMM constructed to date for Pseudomonas. It provides a powerful systems pharmacology platform for the elucidation of complex killing mechanisms of antibiotics.

The exact mode of action of polymyxins is not clear except for the initial electrostatic and hydrophobic interactions with lipid A, a component of the lipopolysaccharide (LPS) in the bacterial outer membrane (OM).Subsequently, the cell envelope is disorganized, cellular contents leak, oxidative stress increases, and finally cell death occurs [2,9,12,13].After polymyxin treatment, P. aeruginosa modifies its lipid A structure to attenuate the aforementioned electrostatic interactions [14].Our recent metabolomics data demonstrated that, apart from lipid A modifications, numerous biochemical pathways are perturbed by polymyxin treatment, indicating that the development of polymyxin resistance by P. aeruginosa involves a complicated interplay of multiple cellular processes [15].There are significant gaps in the knowledge base of the mechanisms of polymyxin activity and bacterial responses in P. aeruginosa, thereby necessitating comprehensive investigations using systems pharmacology approaches.
With the rapid development of genome-scale metabolic models (GSMMs) and the associated flux balance analysis (FBA) methods, systematic investigations into the metabolic changes in response to external nutrient alterations, genetic perturbations, and antibiotic treatments become feasible [16][17][18][19][20][21][22][23][24].Several studies used transcriptomics data as constraints to compute condition-specific metabolic flux changes in response to antibiotic treatments in MDR bacteria, including Acinetobacter baumannii [25], Mycobacterium tuberculosis [26], and Yersinia pestis [27].For P. aeruginosa, 4 GSMMs have been constructed, iMO1056 [28], Opt208964 [29], iMO1086 [30], and the latest iPae1146 [31].iMO1056, Opt208964, and iPae1146 used SEED metabolite and reaction names; iMO1056 and Opt208964 are fully accessible via Model SEED [29,31,32]; iMO1086 used different identifiers (IR/RR plus 5 digits for reactions and C/EC plus 4 digits for metabolites) [30].The previous applications of these models have included simulating the metabolic dynamics in cystic fibrosis patients [33], elucidating the mechanisms of biofilm formation [34,35], predicting potential drug targets [36][37][38], and identifying the key genes that control virulence factors [31].As important as they have been, these models have several overarching limitations.Those past models do not include a major cellular component, the periplasmic space; have poor representation of glycerophospholipid (GPL) biosynthesis; and lack lipid A modification reactions.Considering the pathogenesis of P. aeruginosa, these major limitations significantly compromise the modeling functions.In particular, the power of the 4 reported GSMMs to predict metabolic responses to antibiotic treatment is very limited, as periplasmic GPL and LPS biogenesis plays critical roles in responses to anti-pseudomonal antibiotics such as polymyxins [15,[39][40][41][42].
Here, we describe iPAO1, a newly developed GSMM for P. aeruginosa PAO1 based upon Opt208964 [29] and iMO1056 [28] but with intensive manual curation using several major databases and the literature.Most notably, iPAO1 is the first GSMM for P. aeruginosa where the periplasmic space compartment is incorporated to comprehensively represent cross-membrane transport, GPL metabolism, and LPS biosynthesis.To the best of our knowledge, iPAO1 represents the most comprehensive metabolic reconstruction for Pseudomonas thus far.Modeling with iPAO1 revealed that the lipid A modifications might exert limited impact on cell growth and metabolism but change the physiochemical properties of bacterial OM.Constrained by gene expression levels, the model was used to elucidate the metabolic responses to polymyxin B treatment.Together, iPAO1 provides a powerful systems platform for antimicrobial pharmacological research to combat the rapidly increasing resistance.

Data Description
The genome sequence and annotation of P. aeruginosa PAO1 were obtained from GenBank (accession NC 002516.2).Models iMO1056 and Opt208964 were retrieved from Model SEED [32].The gas chromatography-mass spectrometry (GC-MS) metabolomics data were collected from the literature [43].Metabolites, reactions, and pathways were obtained from databases KEGG (Kyoto Encyclopaedia of Genes and Genomes) [44], MetaCyc [45], TCBD (Transporter Classification Database) [46], TransporterDB [47], and Pseudomonas Genome DB [48].Growth phenotypes on 190 carbon sources and 95 nitrogen sources were determined using BIOLOG Phenotypic Microarrays.Nonessential gene lists were collected from 2 previously reported transposon mutant libraries for PAO1 [49,50].Lipid A of wild-type P. aeruginosa PAK was extracted using a mild acid hydrolysis method, and the structural analysis of lipid A was conducted using mass spectrometry [42].RNA was extracted and used to construct cDNA libraries for RNA-sequencing (RNAseq) on the Illumina MiSeq platform [51].The raw reads were quality trimmed and aligned to the PAO1 reference genome using SubRead [52].Counts were normalized, and the differential gene expression was determined using voom/limma packages with Degust [53].Whole-cell lipids and intracellular metabolites were extracted and analyzed using liquid chromatographymass spectrometry [14,42].Raw lipidomics and metabolomics data were processed with IDEOM software followed by bioinformatic analysis [54].

Strain, media, and BIOLOG experiments
Pseudomonas aeruginosa PAO1 was cultured in LB media and subcultured on nutrient agar.Cells were swapped into sterile capped tubes containing 16 ml IF-0 solution (Cell Biosciences, West Heidelberg, Australia) until the turbidity achieved 42% transmittance in a Turbidimeter (Pacificlab, Blackburn, Australia).The cell suspension was then diluted 5 times with IF-0 solution and dye (Cell Biosciences) to final 85% transmittance.BIOLOG PM 1-3 (Cell Biosciences) were used to investigate the carbon and nitrogen utilization with 2 independent biological replicates.Sodium succinate was used as the carbon source for examining nitrogen utilization.Growth was detected after 24h incubation at 37 • C, using an Infinite M200 microplate reader (Tecan, Mannedorf, Switzerland) at 595 nm.Readings that were ≥1.5-fold of the negative control (i.e., growth media without bacteria) indicated the utilization of nutrients.

Development of a GSMM for P. aeruginosa PAO1
To expedite model development, 2 curated models for PAO1 with the same identifier systems from Model SEED [32], iMO1056 [28], and Opt20896434 [29] were merged.Databases including KEGG [44], MetaCyc [45], and Pseudomonas Genome DB [48] and the literature were used to complete the model with missing components.The identifiers of metabolites and reactions were kept consistent with Model SEED [29] and cross-referred to MetaCyc, KEGG, PubChem [55], ChEBI [56], ChemSpider [57], and BiGG [58].The PAO1 genome annotation from Pseudomonas Genome DB [48] was used to construct "gene to protein to reaction" associations [59].A periplasm compartment was incorporated into the model.Reactions and metabolites were then assigned to cytoplasm, periplasm, and external environment according to the localization prediction of metabolic enzymes by PSORTb 3.0 [60].Transport reactions were generated to enable material exchange across membranes according to TCBD [46] and Trans-porterDB [47].The model was constructed using the Systems Biology Markup Language [61,62].VANTED [63] was used for visualization and analysis of the metabolic network.For each metabolite in the model, specific features, including compartment localization, mass, charge, formula, formation free energy, database identifiers, and source, were added (Additional file 14).Each reaction entered into the model was checked with elementary and charge balance.Reversibility was determined first from the primary literature for each particular enzyme or reaction, if available.Further curation on reaction reversibility and directions was conducted based on change of free energy and knowledge about the physiological direction of a reaction in a pathway.
The Gapfind function from the COBRA toolbox [64] was used to identify the isolated and dead-end metabolites in the model.Candidate reactions from KEGG, MetaCyc, and BiGG were manually inspected for relevance and homology evidence using BLASTp; reactions catalyzed by homologous enzymes (Evalue <1 × 10 −5 , identity ≥35%, coverage ≥50%) were added to the model to eliminate the gaps.Mispredictions of BI-OLOG growth phenotypes were used to refine the draft model (iPAO1 draft2).Further curation was performed to represent the complex biosynthesis pathways of macromolecules (e.g., peptidoglycan, GPL, and LPS).
The biomass formation equation that consisted of the necessary building blocks for bacterial growth was created using the one from iMO1086 [30], with slight modifications on compositions of ions, peptidoglycans, GPL, and LPS (Additional file 17).The growth-and nongrowth-associated maintenance were from iMO1086 [30].

Growth prediction in BIOLOG media
iPAO1 was used to predict the growth phenotypes on chemically defined media with 190 carbon and 95 nitrogen sources (Addi- tional file 18) using the FBA method [24].The objective function of biomass formation was maximized with the specific nutrient uptake rate set at 10 mmol • gDW −1 • h −1 under aerobic condition, as follows: where v biomass denotes the biomass formation flux and S represents the stoichiometric matrix; each metabolic flux v i was constrained by lower and upper bound a i and b i , respectively.All modeling procedures were performed with the COBRA toolbox [64] in MATLAB.The calculated specific growth rates v biomass were then compared to the BIOLOG PM data to assess the prediction accuracy using Fisher's exact test.

Gene essentiality prediction
In silico single-gene deletion was performed using the COBRA toolbox; then, the mutant models were used to predict the specific growth rate in LB broth [32] using FBA.Genes with 99% reduction of the specific growth rate relative to the wild type were defined as essential for cell growth; otherwise, they were considered as semi-essential (1-99% reduction) and nonessential (<1% reduction).Two existing PAO1 transposon insertion mutant libraries, 2-allele mutant library [50,65] and mini-Tn5 insertion mutant library [49], were used to assess the overall prediction accuracy with Fisher's exact test.To determine polymyxin-specific essential genes, transcriptomic constrains were imposed (below) before conducting in silico single-gene deletion simulations.The calculated essential genes identified in polymyxin treatment alone but not in the control were considered as polymyxin specific.

Simulation of bacterial growth and metabolic phenotype changes in response to lipid A modifications
The LPS stoichiometric coefficients in the biomass formula under the control and lipid A modification conditions were set according to the measured lipid A compositions in the wild-type P. aeruginosa PAK in the absence and presence of polymyxin B treatment (Table 1) [14].Aerobic growth was simulated on minimal media with glucose uptake at 10 mmol • gDW −1 • h −1 .For each simulation, the solution space was sampled with 10000 random points using the ll-ACHRB algorithm [66].Flux samples of the control and lipid A modification were then compared.Significantly perturbed metabolic fluxes were identified using a Zscore based approach [67].
To further analyze the metabolic impact of lipid A modifications, the proportions of all types of LPS in the biomass formula were randomly assigned and the process was repeated 1000 times.For each repetition, the specific growth rates were calculated and solution space was sampled using the methods described above.For each type of lipid A, specific physiochemical properties (f) including total atom number, partition coefficient (logP), average charge, and molecular polarity were predicted at pH 7 using the cxcalc tool from ChemAxon (Budapest, Hungary).The overall apparent properties F of the OM were estimated by calculating the weighted sum, as follows: where w j represents the stoichiometric coefficient of the j-th of 288 heterogeneous LPS molecules in the biomass formula.Pairwise correlation analysis was conducted between lipid A modifications, physiochemical properties changes, bacterial growth, and metabolism alterations.

Predict metabolic responses to polymyxin treatment by constraining fluxes with transcriptomics data
The RNA-seq data from 1-h 1 mg•L −1 polymyxin B treatment experiment using PAO1 were used as flux constraints for modeling [51].For each gene under every condition, the reads per kilobase million (RPKM) value was calculated from the aligned reads using the edgeR package [68] and normalized to constrain flux upper bounds (b i ) using the E-Flux algorithm [26].Specifically, for each reaction catalyzed by a single enzyme, the upper flux bound was set to the determined RPKM value under the respective condition.For a reaction catalyzed by an enzyme complex, the upper bound was set to the minimum RPKM value of the associated genes.For a reaction catalyzed by isozymes, the upper bound was set to the sum of RPKM values of the associated genes.The maximum of upper bounds was then normalized to 10000 mmol • gDW −1 • h −1 .The lower bounds a i were set to 0 for irreversible and −b i mmol • gDW -1 • h −1 for reversible reactions, respectively.CAMHB was used in the RNA-seq experiment; it is known as an undefined medium that contains mainly amino acids and oligopeptides [69].The maximum uptake rates of amino acids in P. aeruginosa vary between 0.26 and 1.44 mmol [70,71,73].Therefore, the upper bounds (b CAMHB i ) of uptake rates of amino acids, vitamins, and dipeptides in iPAO1 were constrained to 1 mmol • gDW −1 • h −1 without loss of generality.For each condition, the solution space was sampled with 10000 points using ll-ACHRB as described above.Statistical significance of differential flux distributions was computed using the Z-score method described above.The turnover rate for key metabolites was calculated by summing up all influxes or effluxes [74].To assess the impact of changing nutrient uptake bounds, sensitivity analysis was conducted by randomly sampling solution space as above while varying b CAMHB i from 0.26 to 1.44 mmol
Second, a periplasmic compartment was built to incorporate 698 periplasmic metabolites, 509 transport reactions across the inner membrane (IM), 441 transport reactions across the outer membrane (OM), and 403 periplasmic reactions.The resulting intermediate model was designated as iPAO1 draft2.
Third, the major pathway gaps were filled.GapFind [75] identified 109 dead-end metabolites (Additional file 7).The growth phenotypes on 190 carbon and 95 nitrogen nutrients were predicted using iPAO1 draft2, and compared with our experimental BIOLOG phenotypic microarray (PM) results (Additional file 8).As a result, 162 false-negative predictions (i.e., the prediction indicated nongrowth whereas the BIOLOG experiment demonstrated valid growth on a specific nutrient) were determined, indicating the lack of associated transport or catabolic reactions for these nutrients.To link the dead-end metabolites back to the metabolic network and eliminate inconsistencies with the BIOLOG PM results, several modifications were made, including adjustment of the reversibility settings of 180 reactions and change in the directions of 87 reactions (Additional file 9); removal of 14 metabolites and 96 reactions (Additional files 10 and 11), which were either duplicated (e.g., β-D-glucose was duplicated with D-glucose) or representing general metabolite classes (e.g., protein, mRNA, DNA); and addition of 98 boundary reactions, 677 transport reactions, and 252 metabolic reactions (Additional file 12).Resolving the false-negative predictions of the BIOLOG growth phenotypes substantially improved the model.For example, predictions using iPAO1 draft2 showed that PAO1 was unable to grow with formic acid as a sole carbon source due to lack of the corresponding transport reaction.Interrogation of the Pseudomonas Genome Database [48] and Pfam [76] identified PA2777, a hypothetical protein in the National Center for Biotechnology Information (NCBI) and UniProt that may encode formic/nitrite transporter (Pfam01226, P = 7e-34).Subsequent addition of the transport reaction (rxn08526) enabled in silico growth of PAO1 on formic acid.Another example is that initially iPAO1 draft2 failed to predict utilization of 1,2-propanediol for growth owing to the exiting gap in dehydrogenation of 1,2-propanediol to lactaldehyde.Using the Basic Local Alignment Search Tool for Proteins (BLASTp) with the query sequence of lactaldehyde reductase (fucO, b2799) from Escherichia coli K12 MG1655, we identified a candidate homologous gene PA1991 (Identity = 35%, Eval = 2e-75, BLASTp).PA1991 encodes an iron-containing alcohol dehydrogenase and has more than 300 orthologues in Gram-negative bacteria that encode lactaldehyde oxidoreductases or 1,2-propanediol dehydrogenases according to OrthoDB [77].Inactivation of PA1991 resulted in an 8-fold prolonged lag phase when P. aeruginosa grew on 1,2-propanediol [78].Therefore, reaction rxn01615 oxidizing 1,2-propanediol to lactaldehyde was added into iPAO1 draft2.A very large number of such labor-intensive manual curations were conducted to improve the model.This enabled in silico growth on a number of nutrients from BIOLOG PMs, including 4hydroxyphenylacetate, tyramine, quinic acid, itaconic acid, citramalic acid, L-pyroglutamic acid, carnidine, glycinebetaine, Lmethylsuccinate, and D-amino acids (Additional file 8).
Fourth, the biogenesis of the bacterial envelope was delineated.Cross-linking between amino acids residues among peptidoglycan chains results in a rigid network structure in P. aeruginosa [79].In total, 17 reactions representing peptidoglycan cross-linking and hydrolysis were incorporated by searching for homologues of glycosyltransferases, transpeptidases, carboxypeptidases, and endopeptidases in PAO1 [80].Overall, a

Inner membrane
Outer membrane detailed peptidoglycan biosynthesis pathway was constructed with 60 reactions.GPL compositions in the bacterial membranes can change in response to antibiotic treatment [39,81].Previous studies [82] and our own lipidomics results [14] showed a great diversity in GPL species in P. aeruginosa.Overall, 386 unique metabolites (i.e., 66.2% of the 583 metabolites in the GPL metabolism pathway) and 367 reactions (66.7% of the 550 reactions in the GPL metabolism pathway) were incorporated into iPAO1 draft2 (Additional files 1, 13, and 14; Fig. 1).LPS consists of lipid A, core oligosaccharide, and O-antigen polysaccharide [40] and plays key roles in the host-pathogen interaction and resistance to antibiotics such as polymyxins [13,83].A detailed synthesis and interconversion network was generated with 432 types of LPS and 1169 reactions (Fig. 2; Additional file 1).Notably, our GSMM is the most comprehensive to date in lipid A biosynthesis and modifications.
Components in iPAO1 were aligned with databases including KEGG [44], MetaCyc [45], PubChem [55], ChemSpider [57], ChEBI [56], Model SEED [32], and BiGG [58] (Additional files 15 and 16).Consequently, 1404 (46.5%), 1590 (52.6%), and 2142 (70.9%) metabolites have corresponding identifiers in MetaCyc, KEGG, and Model SEED, respectively; 1556 (35.6%), 1596 (36.6%), and 1964 (45.0%) reactions were computationally mapped to the reactions from MetaCyc, KEGG, and Model SEED, accordingly (Fig. 3B).A significant portion of mismatches were caused by the incorporation of specific types of metabolites in the GPL metabolism and LPS biosynthesis pathway, which in databases are usually lumped as general compound classes.The properties of metabolites, including mass, charge, and formula, were included in iPAO1.The standard Gibbs free energy changes of formation ( f G • ) and reaction ( r G • ) were obtained from  A breakdown of genes involved in iPAO1 (Additional file 17) using the clusters of orthologous groups (COGs) showed remarkable improvement compared to previous reconstructions (Fig. 3C).The largest increase in the coverage compared to iMO1056 is lipid transport and metabolism (24.1%), followed by inorganic ion transport and metabolism (19.3%); whereas compared to Opt208964, the largest increase in the coverage is nucleotide transport and metabolism (57.9%), followed by amino acid transport and metabolism (52.0%).Overall, the transport and metabolism of nucleotides and amino acids showed the highest percent coverage of COG functional categories in iPAO1 (72.9% and 65.6%, respectively).Notably, the reactions in categories not apparently related to metabolism were dramatically reduced in iPAO1 compared to Opt208964, including translation, ribosomal structure and biogenesis, post-translational modification, protein turnover, chaperones and signal transduction mechanisms, and undetermined categories, including function unknown class.
In iPAO1, GPL metabolism, LPS biosynthesis, and transport across OM were ranked the 3 largest pathways and also contained the largest proportion of curated reactions (Fig. 3D).Additionally, these 3 pathways have high reaction-to-gene ratios (13.1-24.2;Fig. 3E), indicating that enzymes in these pathways are capable of acting on a broad range of substrates.As kinetic parameters are usually not involved in a GSMM, constraintbased analyses (e.g., FBA) of a GSMM do not directly account for enzyme levels, intracellular metabolic concentrations, or substrate-level regulation.Accordingly, the affinity difference of various substrates was not considered in our iPAO1 modeling effort.
We used the biomass formation equation from iMO1086 to construct iPAO1 with modifications on LPS and ion species (Additional file 18).In addition, to take into account the ex-tra energy consumption caused by charging tRNAs, the original amino acids in the biomass formation reaction were replaced by aminoacyl-tRNA, followed by addition of specific charging reactions to the model.Taken together, iPAO1 represents the most comprehensive metabolic reconstruction thus far for P. aeruginosa PAO1.

Growth capability on various nutrients
Investigation of nutrient utilization using BIOLOG PMs showed that PAO1 could utilize a broad range of nutrient sources, indicated by the observed growth on 68 of 190 (35.8%) carbon and 76 of 95 (80.0%) nitrogen substrates (Fig. 4).Growth simulation with iPAO1 achieved an overall accuracy of 89.1% (254 of 285), which substantially outperformed previous models (81.5% for Opt208964 [29], 77.9% for iMO1056 and iMO1086 [30], and 80% for iPae1146 [31]).Twenty-one false-positive and 10 false-negative (Fig. 4, Additional file 8) disagreements were observed, possibly due to the complexity of regulatory mechanisms and missing annotation of nutrient transport and/or catabolism pathways in PAO1.

Impact of lipid A modifications on bacterial growth and metabolism
Pseudomonas aeruginosa modifies lipid A components in the OM in response to polymyxin treatment [85].The LPS stoichiometric coefficients in the biomass formula of iPAO1 were  No. of components ; the generated (R)-3-hydroxydecanoate was fuelled into β-oxidation to produce octanoyl-CoA, which was subsequently salvaged for fatty acid biosynthesis.
To further investigate the impact of lipid A modifications on bacterial growth, 1000 sets of the compositions of 288 hetero-geneous LPS molecules were randomly generated with the total proportion of LPS unchanged in the biomass formation formula (Additional file 21).The metabolic fluxes were calculated for each of the 1000 sets of LPS compositions by randomly sampling the solution space with 10000 points.Across the 1000 sets of metabolic fluxes (Additional file 23), the specific growth rate varied between 0.8812 and 0.8897 mmol • gDW −1 • h −1 .Correlative analysis of the apparent overall physiochemical properties of lipid A (Additional file 22) with the predicted growth phenotypes showed 3 interesting findings.First, addition of L-Ara4N reduced the negative charge of lipid A (ρ = 1.00), decreased the hydrophobicity of the OM (represented by logP, ρ = −0.59)but required assimilation of more ammonia (represented by ammonia turnover, ρ = 0.57).Second, hydroxylation on acyl chains of lipid A exerted minor effects over either bacterial growth or physiochemical properties.Third, addition of acyl chains resulted in large lipid A molecules (represented by the atomic counts,  ρ = 0.88), enhanced molecular polarity of lipid A (ρ = 0.87), increased OM hydrophobicity (ρ = 0.75), and, notably, retarded growth (ρ = −0.95),reduced redox and energy turnover (ρ = −0.98 for both), and increased requirement of ammonia (ρ = 0.59) (Fig. 5).It is evident that none of the three aforementioned modifications produced a dramatic impact on bacterial growth or metabolism (Additional file 23).

Elucidating the mechanisms of metabolic responses to polymyxin treatment
RNA-seq data were utilized as model constraints (Additional file 24) with an E-Flux method [67] to calculate the metabolic fluxes in the absence and presence of polymyxin B (see the Methods section).The exchange fluxes were constrained based on the maximum uptake rates of the media ingredients (see the Methods section and Additional file 1).Comparison of the flux distributions revealed that 1392 reactions were differentially regulated (FDR <0.01,Additional file 25).A range of metabolic pathways were significantly disturbed, including central metabolism, amino acid metabolism, purine biosynthesis, fatty acid biosynthesis and metabolism, LPS and GPL biosynthesis, and transport reactions.Polymyxin B treatment reduced the growth rate (18.2%) and increased oxygen uptake (6.9%) and CO 2 emission (6.0%); however, the respiration quotient remained roughly unchanged (Table 3).
Calculating the flux-sum of critical cofactors revealed 13.1% increase of redox turnover and 8.2% decline of energy turnover Table 3: Specific growth rate, significantly altered major exchange fluxes (>1 mmol • gDW −1 • h -1 ), respiration quotient, and the fluxes through F 0 F 1 -ATPase calculated using the RNA-seq data [51]  a FDR was calculated using the Benjamini-Hochberg method [72].
after 1 mg • L −1 polymyxin B treatment for 1 h.Breaking down the cofactors showed that the turnover of major redox equivalents NADH, Reduced nicotinamide adenine dinucleotide phosphate (NADPH), ubiquinol-8, and FADH2 substantially increased by 12.6%, 13.9%, 3.9%, and 35.9%, respectively; whereas the turnover of ATP, the major contributor to energy, significantly decreased by 8.52% after 1 mg • L −1 polymyxin treatment for 1 h (Fig. 6, Additional file 26).Overall, metabolic flux analysis using iPAO1 integrated with our transcriptomics data revealed a significant global impact on bacterial metabolism due to polymyxin B treatment.

Discussion
The emergence of Gram-negative 'superbugs' that are resistant to the last-resort polymyxins highlights the urgent need for novel approaches such as GSMMs to understand the mechanisms of antibacterial activity and resistance.The main utility of GSMMs is their ability to bridge critical gaps between genomics and metabolic phenotypes through the prediction of metabolic responses to antimicrobial treatments at the network level.Here, we report the development, optimization, validation, and application of a high-quality GSMM designated iPAO1 for a type strain P. aeruginosa PAO1.Importantly, iPAO1 was used to understand the complicated effect of polymyxin treatment on bacterial metabolism.Simulation with iPAO1 showed that lipid A modifications in response to polymyxin treatment only exert minor effects on bacterial growth and metabolism.Further calculations that integrate transcriptomics data as model constraints revealed that polymyxin treatment may reduce growth and affect a broad range of pathways.iPAO1 represents the most comprehensive metabolic model for P. aeruginosa to date and incorporates 1458 genes, accounting for approximately 25.8% of the PAO1 genome.Among the 4 GSMMs developed for P. aeruginosa PAO1, iMO1086 and iPae1146 were constructed on the basis of iMO1056 with moderate increase of metabolites, reactions, and genes [28,30,31]; Opt208964 is also in a medium size, which limits modeling capacity [29].In contrast, iPAO1 is significantly expanded in model scale by doubling or even tripling the numbers of metabolites and reactions (Fig. 3A).iPAO1 achieved an unprecedented prediction accuracy of 89.1% for growth on various nutrients, out-performing all of the previously reported GSMMs for P. aeruginosa [28][29][30][31].The iPAO1 model was also used to predict gene essentiality with a high accuracy of 87.9%.Given the extensive curation and significant expansion, iPAO1 will serve as the primary reference for future development of metabolic models, particularly for other P. aeruginosa strains.
Unlike iPAO1, none of the previous P. aeruginosa GSMMs incorporated the periplasm.As polymyxins initially target LPS in the OM and can cause substantial changes in the cell envelope, the periplasmic space is a major component in iPAO1.The periplasmic space of E. coli is estimated to constitute up to 16% of total cell volume [86].It contains a thin cell wall composed of peptidoglycan and a variety of ions and proteins, which are involved in transport, folding, cell envelope biogenesis, electron transport, and xenobiotic metabolism [87].iPAO1 is the first P. aeruginosa GSMM to incorporate the periplasmic compartment, enabling accurate representation of metabolic machinery, especially for those reactions that occur exclusively in this important cellular space, and transport of substrates across the IM and OM.Furthermore, iPAO1 provides detailed representations of GPL and LPS biosynthesis that allows the precise mapping of GPL and LPS responses from experimental metabolomics and lipidomics data (Fig. 1 and 2).
In response to polymyxin treatment, Gram-negative bacteria modify their lipid A with cationic moieties (i.e., phosphoethanolamine and L-Ara4N) that act to repel the like-charge of the polymyxin molecule [40].Based on our simulations (Additional file 20), we purport that such lipid A modifications exerted a limited impact on cellular metabolism and growth.Most of the flux changes were insignificant; the remaining significant flux changes mainly resulted from futile cycles that contain sets of reactions using redox equivalents, whereas the net carbon flow remained unchanged.Simulation using randomized lipid A compositions further consolidated our hypothesis that lipid A modifications cause only moderate variations of bacterial growth and metabolism (Fig. 5, Additional file 23).Our simulation results revealed that lipid A modifications result in substantial physiochemical changes in the OM of P. aeruginosa, including neutralizing the surface negative charge by addition of positively changed L-Ara4N and altering the polarity and hydrophobicity by acylation and deacylation.The general mode of action of polymyxin involves the initial electrostatic interaction between the cationic side chains of the polymyxin molecule with the anionic lipid A head groups [83].These events are subsequently followed by hydrophobic interactions between the N-terminal fatty acyl chain and position 6/7 hydrophobic side chains of the polymyxin with the hydrophobic fatty acyls of lipid A [83].Therefore, in concept, both the addition of L-Ara4N and deacylation of lipid A should contribute to polymyxin resistance.Indeed, in our recent transcriptomic and neutron reflectometry studies, we discovered that deletion of the corresponding gene pagL (PA4661) resulted in an increased susceptibility to polymyxins in a polymyxin-resistant mutant PAKpmrB6 derived from P. aeruginosa PAK [14,88], demonstrating that the lipid A deacylation also plays a key role in the response of P. aeruginosa to polymyxin treatment.Also, in our recent transcriptomics and metabolomics studies, we discovered that polymyxin treatment leads to remarkable growth reduction and metabolic perturbations in Gram-negative bacteria [41,42,[89][90][91].The integration of transcriptomics results into GSMMs allows for more accurate predictions of metabolic responses to either environmental (i.e., antibiotic treatment) or genetic perturbations (i.e., mutations) [92].In the present study, we used the E-Flux method to integrate transcriptomics data as flux constraints [26].E-Flux can map continuous gene expression levels to the metabolic network and uses the transcript abundance to determine the degree to which a reaction is active or inactive [26].Therefore, E-Flux provides a more Zhu et al.
physiologically relevant description of the continuous nature of the reaction activity and avoids use of any artificial thresholds to binarize gene expression data [93].In the present study, metabolic fluxes with and without antibiotic treatment were not calculated with minimization of metabolic adjustment (MOMA), as MOMA was developed to predict the metabolic flux redistributions in gene knockout mutants [94].MOMA hypothesizes that metabolism of the mutant tends to approximate the wild type [94], which is distinct from the antibiotic treatment scenario.For instance, our metabolomics data have demonstrated that the antibiotic treatment caused dramatic metabolic changes in bacteria [41].
Comparison of the calculated flux distributions revealed that a broad range of metabolic perturbations occur in response to polymyxin treatment (Fig. 6), ranging from central carbon metabolism to oxidative phosphorylation and amino acid metabolism.Reduced growth, increased redox turnover, and decreased energy turnover due to polymyxin treatment were evident (Fig. 6), indicating that bacterial cells regulated their metabolism to produce more redox power to cope with the oxidative stress.This is consistent with previous findings that showed bactericidal antibiotics induced lethal oxidative damages by generating highly deleterious free radicals with subsequent culmination of cellular death [95].In addition, our simulations revealed that polymyxin treatment induced an uptake of L-alanine, which was catabolized to generate more NADH (Fig. 6).This indicates that rich media (e.g., CAMHB) may provide abundant amino acids and peptides that can be utilized by bacterial cells to generate sufficient redox equivalents to cope with the oxidative damage caused by polymyxin treatment.Furthermore, our simulation results also showed an upregulated metabolic flux toward L-spermidine biosynthesis upon polymyxin B treatment (rxn00127 and rxn01406; Additional file 25).Previous studies showed that polyamines (e.g., spermidine) could protect P. aeruginosa from antimicrobial peptide killing [96].It is assumed that the cationic spermidine could interact with the anionic LPS, mask the negative cell surface, and reduce the electrostatic interactions between polymyxin B and bacterial Therefore, the enhanced biosynthesis of spermidine might increase its abundance at the cell surface and contribute to polymyxin resistance.
The constructed iPAO1 provides a detailed presentation of LPS biogenesis (Fig. 2), in particular, lipid A modifications.Further integration with specific regulatory modules will enable dynamic simulation of metabolic responses to polymyxin treatment.Previous studies revealed that various two-component regulatory systems (TCSs), including PhoPQ, PmrAB, ParRS, CprRS, and ColRS, play key roles in regulating polymyxin resistance [84,[97][98][99][100].Among them, the PmrAB and PhoPQ systems are able to sense the depletion of external cations (e.g., Mg 2+ and Ca 2+ ) and upregulate the expression of the arnBCADTEFugd operon, which is responsible for the modification of lipid A with L-Ara4N [101].Moreover, the fatty acylation of lipid A by PagP is under the control of PhoPQ [102,103].ParRS and CprRS are independent TCSs that mediate the upregulation of pmrAB, arnBCADTEF-ugd operon, pagL, and adaptive resistance in response to polymyxin treatment [97,104].In overview, lipid A modifications due to polymyxin treatment are strictly controlled by very complex regulatory networks that involve signal sensors, transcriptional regulators, and metabolic enzymes.Therefore, future studies are warranted to integrate these regulatory modules into the GSMM to enable simulating bacterial response dynamics to polymyxin treatment and analyzing adaptive resistance mechanisms in P. aeruginosa.
Overall, we have constructed, optimized, and validated a high-quality genome-scale metabolic model iPAO1 for P. aeruginosa PAO1.This comprehensive model incorporates metabolic pathways, particularly the biogenesis of membrane components, and enables delineation of the complex metabolic responses to antibiotics.iPAO1 provides a valuable systems tool for quantitative simulation of bacterial metabolic responses to antibiotics, elucidation of the molecular mechanisms of antimicrobial killing and resistance, and facilitation of designing rational antimicrobial combination therapy.To the best of our knowledge, this study is the first to integrate antimicrobial pharmacology, computational biology, a metabolic network, and systems pharmacology to analyze large-scale datasets in order to better understand the dynamic and complex nature of polymyxin killing and resistance.Combined with antibiotic pharmacokinetics and pharmacodynamics, iPAO1 offers an in silico platform for precision polymyxin chemotherapy.

Conclusion
The generated collection of transcriptomics, metabolomics, lipidomics, and lipid A profiling data provides comprehensive datasets of P. aeruginosa for future integrative analysis of polymyxin systems pharmacology.As the largest curated GSMM thus far for Pseudomonas, iPAO1 represents all aspects of the cellular metabolism and may serve as the platform for integrative analysis of multiomics data.Simulation with transcriptomics constraints in this study revealed metabolic flux changes in amino acid catabolism, tricarboxylic acid cycle, and redox turnover caused by polymyxin treatment.Correlative analysis of metabolomics and transcriptomics with the constraint-based modeling is necessary for delineating the regulatory effects on metabolism.The methodology of using GSMMs to analyze multilevel omics data is applicable to other areas beyond antimicrobial pharmacology.Further integration with antimicrobial pharmacokinetics and pharmacodynamics will not only provide better pharmacological understanding but also empower the model to quantitatively predict the bacterial responses to antimicrobial therapy in the context of complex interplay of signaling, transcriptional regulation, and metabolism.In summary, our GSMM approach provides a powerful systems tool to elucidate the complex mode of action of antibiotics and will paradigm shift antimicrobial pharmacology.

Availability of supporting data
The raw RNA-seq data have been deposited in the NCBI sequence read archive database under the BioProject accession number PRJNA414673.The metabolomics and lipidomics data have been deposited in the Metabolight database with accession number MTBLS630.Supporting data, including the scripts used in this project, are available via the GigaScience repository GigaDB [105].Additional file 1 (additionalFile1.docx):Manual curation of GPL biosynthesis, LPS biosynthesis and modification pathways, and sensitivity analysis of nutrient uptake bounds.

Figure 2 :
Figure 2: LPS biosynthesis and modification in iPAO1.A) VANTED diagram showing the biosynthesis of different LPS molecules.B) LPS biosynthesis pathway; lipid A and LPS are indicated in the same color as in (A).

Figure 3 :
Figure 3: Constitutional genes, reactions, and metabolites in iPAO1.A) Sources of iPAO1 components.B) Radar map showing the percentages of metabolites and reactions with valid database identifiers.C) The COG functional classification of the involved genes in iMO1056, Opt208964, and iPAO1.Percentages given in the middle indicate the coverages of COG groups.The proportions of the curated reactions (D), reaction-to-gene ratio (E), and predicted subcellular localizations of the involved proteins (F) are shown for each pathway or COG group.In panel D, red bars indicate the curated reactions; blue bars indicate the reactions from the previous model.In panels D and E, pathways with the highest curation proportion or reaction-to-gene ratio are highlighted in red.

Figure 4 :
Figure 4: Comparison of the BIOLOG result (left columns) and model prediction (right columns).Blue indicates growth and yellow indicates no growth.

Figure 5 :
Figure 5: Simulation of the impact of lipid A modifications on bacterial growth, metabolism, and OM physiochemical properties.The significant correlation (P < 0.05) of paired items is indicated in red.

Figure 6 :
Figure 6: Polymyxin B-induced metabolic perturbations.The distributions of metabolic fluxes and metabolite turnover rates are shown in subgraphs with red indicating control and blue indicating polymyxin B treatment.

Figure S1 :
Figure S1: Sensitivity analysis of the mean metabolic fluxes (A) and metabolite turnover rates to the variation of nutrient uptake upper bounds.Red indicates the control and blue indicates polymyxin B treatment.Additional file 1 (additionalFile1.docx):Manual curation of GPL biosynthesis, LPS biosynthesis and modification pathways, and sensitivity analysis of nutrient uptake bounds.

Table 1 :
[14]d A composition (%) in the outer leaflet of the OM in PAK[14]

Table 2 :
Components in model iMO1056, Opt208964, and iPAO1 a Pathway information is not available in iMO1056 from the Model SEED database.MetaCyc and Model SEED for 1877 metabolites (62.1%) and 1355 reactions (31.0%) (Additional files 15 and 16).