Next Article in Journal
Longitudinal Characterization of the Gut Microbiota in the Diabetic ZDSD Rat Model and Therapeutic Potential of Oligofructose
Previous Article in Journal
Computational Insights into Natural Antischistosomal Metabolites as SmHDAC8 Inhibitors: Molecular Docking, ADMET Profiling, and Molecular Dynamics Simulation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Holistic Approach from Systems Biology Reveals the Direct Influence of the Quorum-Sensing Phenomenon on Pseudomonas aeruginosa Metabolism to Pyoverdine Biosynthesis

by
Diana Carolina Clavijo-Buriticá
1,*,
Catalina Arévalo-Ferro
1 and
Andrés Fernando González Barrios
2
1
Grupo de Comunicación y Comunidades Bacterianas, Departamento de Biología, Universidad Nacional de Colombia, Carrera 45 No. 26-85, Bogotá 111321, Colombia
2
Grupo de Diseño de Productos y Procesos (GDPP), Departamento de Ingeniería Química y de Alimentos, Universidad de los Andes, Edificio Mario Laserna, Carrera 1 Este No. 19ª-40, Bogotá 111711, Colombia
*
Author to whom correspondence should be addressed.
Metabolites 2023, 13(5), 659; https://doi.org/10.3390/metabo13050659
Submission received: 21 March 2023 / Revised: 26 April 2023 / Accepted: 6 May 2023 / Published: 16 May 2023
(This article belongs to the Section Bioinformatics and Data Analysis)

Abstract

:
Computational modeling and simulation of biological systems have become valuable tools for understanding and predicting cellular performance and phenotype generation. This work aimed to construct, model, and dynamically simulate the virulence factor pyoverdine (PVD) biosynthesis in Pseudomonas aeruginosa through a systemic approach, considering that the metabolic pathway of PVD synthesis is regulated by the quorum-sensing (QS) phenomenon. The methodology comprised three main stages: (i) Construction, modeling, and validation of the QS gene regulatory network that controls PVD synthesis in P. aeruginosa strain PAO1; (ii) construction, curating, and modeling of the metabolic network of P. aeruginosa using the flux balance analysis (FBA) approach; (iii) integration and modeling of these two networks into an integrative model using the dynamic flux balance analysis (DFBA) approximation, followed, finally, by an in vitro validation of the integrated model for PVD synthesis in P. aeruginosa as a function of QS signaling. The QS gene network, constructed using the standard System Biology Markup Language, comprised 114 chemical species and 103 reactions and was modeled as a deterministic system following the kinetic based on mass action law. This model showed that the higher the bacterial growth, the higher the extracellular concentration of QS signal molecules, thus emulating the natural behavior of P. aeruginosa PAO1. The P. aeruginosa metabolic network model was constructed based on the iMO1056 model, the P. aeruginosa PAO1 strain genomic annotation, and the metabolic pathway of PVD synthesis. The metabolic network model included the PVD synthesis, transport, exchange reactions, and the QS signal molecules. This metabolic network model was curated and then modeled under the FBA approximation, using biomass maximization as the objective function (optimization problem, a term borrowed from the engineering field). Next, chemical reactions shared by both network models were chosen to combine them into an integrative model. To this end, the fluxes of these reactions, obtained from the QS network model, were fixed in the metabolic network model as constraints of the optimization problem using the DFBA approximation. Finally, simulations of the integrative model (CCBM1146, comprising 1123 reactions and 880 metabolites) were run using the DFBA approximation to get (i) the flux profile for each reaction, (ii) the bacterial growth profile, (iii) the biomass profile, and (iv) the concentration profiles of metabolites of interest such as glucose, PVD, and QS signal molecules. The CCBM1146 model showed that the QS phenomenon directly influences the P. aeruginosa metabolism to PVD biosynthesis as a function of the change in QS signal intensity. The CCBM1146 model made it possible to characterize and explain the complex and emergent behavior generated by the interactions between the two networks, which would have been impossible to do by studying each system’s individual components or scales separately. This work is the first in silico report of an integrative model comprising the QS gene regulatory network and the metabolic network of P. aeruginosa.

1. Introduction

In recent years, computational modeling of biological systems has acquired a significant role in understanding and predicting cellular performance. Computational modeling for explaining cell behavior has been accomplished by numerous discrete, continuous, and spatially explicit models describing subcellular systems and processes that biological systems need to survive [1]. Indeed, all biological systems are complex and dynamic and can be modeled as clusters of various functional networks. Metabolic, cell signaling, gene regulation, and intercellular communication networks are interconnected, regulate each other, and operate in diverse temporal and spatial domains to maintain living organisms’ growth, development, and reproductive potential [2]. The functional networks of biological systems can operate on characteristic and specific temporal and spatial scales. For instance, intracellular biochemical reactions can ensue on time scales of seconds or less, involving gradients exceeding microns in length. Functional phenotypes at the cellular and cell-to-cell levels unfold in tens of seconds to minutes and on spatial scales of tens of microns [3]. Therefore, biological systems are inherently considered as hierarchical and multiscale in time and space. Most traditional computational models focus on studying biological systems at single biological scales. However, systemic approaches facilitate the designing of multiscale computational models, to bridge those scales, so that the effects of perturbations at one scale can be predicted at other scales. Thus, multiscale computational models are exceptionally positioned to capture the relationship between different scales of biological function; they can bridge the gap between isolated in vitro experiments and whole-organism in vivo models [1,2,3].
Following the development of next-generation sequencing techniques, a key objective has been to relate annotated genome sequences to cellular physiological functions. Therefore, systems biology has been designing new computational strategies to reconstruct biological networks (genomic and metabolic) and to model and dynamically simulate biological systems in order to study, from a holistic perspective, the regulation of mechanisms leading to the expression of diverse phenotypes. Thus, the extensive dataset of measurements of all biomolecules at the system level favors integration from the molecular to the whole organism level. In recent decades, several computational studies have made developing multiscale models of some biological systems possible. For example, Thiele et al. and Biggs and Papin have shown the advantages of multiscale models for understanding metabolic and macromolecular synthesis in Escherichia coli and biofilm formation by Pseudomonas aeruginosa, respectively [4,5].
Following a holistic approach, this work aimed to explain the synthesis of the siderophore pyoverdine (PVD), a virulence factor in P. aeruginosa [6]. The PVD synthesis requires two biological processes. The first is the expression of pvd genes, which encode enzymes involved in PVD synthesis and which are regulated by the quorum-sensing (QS) phenomenon [7] (a mechanism of cell communication mediated by signaling molecules -in this case, induced under conditions of iron deficiency, and responsible for synchronizing phenotype expression in a bacterial community—in this case, the factor virulence PVD expression) [8]. The second process corresponds to the metabolic pathway involved in PVD synthesis. Therefore, to understand how these two biological processes work and how they regulate each other, it was necessary to design a methodology to develop a model capable of describing their interactions.
P. aeruginosa has two QS systems, LasI/LasR and RhlI/RhIR. The signaling molecule of the LasI/LasR system is N-(3-Oxododecanoyl)-L-homoserine lactone (3O-C12-HSL), which is synthesized by LasI and detected by the LasR receptor. The RhlI/RhIR system is regulated by N-butanoyl-L-homoserine lactone (C4-HSL), synthesized by RhII, and detected by the RhIR receptor. Both QS systems are hierarchically organized, such that the LasI/LasR system regulates the function of the RhlI/RhIR system, allowing further fine-tuning of the QS feedback on specific genes [9]. These QS systems are linked through regulatory connections involving regulatory proteins (RsaL, PqsR, RclA, vfR, and GacA) and the Pseudomonas Quinolone Signal (PQS) system. This system encodes the synthesis of 2-heptyl-3-hydroxy-4-quinolone (2C7-3OH-4(1H)Q) that forms a complex with the available ferric ion (Fe3+) and plays an essential role in the pvd genes expression. At the cellular level, the quinolone complex positively controls the expression of the Pvd-type proteins, which is responsible for PVD maturation, synthesis, and exportation [10].
Most PVD biosynthesis is carried out in the bacterial cytosol by four non-ribosomal peptide synthase (NRPS) proteins encoded by the pvdL, pvdI, pvdJ, and pvdD genes [6]. Together, these genes encode ferribactin synthase, responsible for synthesizing ferribactin, the precursor of PVD. PVD biosynthesis also requires the expression of the pvdA, pvdF, and pvdH genes. According to the model proposed by Visca et al. [8], the PvdE protein (an ABC-type transporter), located in the inner membrane, transports ferribactin to the periplasmic space where PVD chromophore maturation occurs by an oxidation-reduction mechanism driven by the PvdM, PvdP, PvdQ, PvdO, and PvdN proteins [7,11,12] through a mechanism not fully elucidated. Finally, mature PVD is secreted into the extracellular space by an unclear process that chelates iron and forms the PVD-Fe3+ complex, which is transported back into the intracellular space by the FpvA carrier. Iron is then released from PVD to the periplasm by a mechanism involving iron reduction (Fe3+ to Fe2+), and PVD is then transported out of the cell via a specialized carrier of the siderophore recycling process [13,14,15,16].
An integrative (multiscale) model that combines the QS gene regulatory network and the metabolic network of P. aeruginosa is proposed here for the first time to describe and understand the influence of the QS phenomenon on PVD biosynthesis in P. aeruginosa from a systemic perspective. That model involves the process of bacterial communication mediated by the QS phenomenon for synthesizing signal molecules that regulate the expression of pvd genes to generate enzymes that metabolically catalyze the production of the virulence factor. Furthermore, this integrative model combined two biological network models: (i) The QS gene regulatory network model, based on deterministic approaches and whose simulation results were the concentration of chemical species and the fluxes of chemical reactions; and (ii) the model of the P. aeruginosa metabolic network based on the flux balance analysis (FBA) approximation. This approximation addressed the biomass maximization problem constrained by mass balance and thermodynamic conditions to obtain the optimal distributions of reaction fluxes through the metabolic network [17,18,19,20,21]. The methodological strategy first identified reactions shared by the QS and the P. aeruginosa metabolic networks. Then, the fluxes of these shared reactions, obtained from the QS network, were used as a constraint system to model the P. aeruginosa metabolic network by dynamic flux balance analysis (DFBA) [22,23]. The proposed multiscale model showed a capability to infer the influence of QS on the P. aeruginosa strain PAO1 metabolism to PVD biosynthesis as a function of variations in QS signal intensity.

2. Materials and Methods

Computational modeling followed three steps: (i) construction of a deterministic model of the P. aeruginosa QS network as a gene regulatory network for the expression of pvd genes encoding enzymes for PVD synthesis and maturation; (ii) construction of the P. aeruginosa genome-scale metabolic network, including PVD biosynthetic reactions using the FBA approximation; and (iii) integration of these two network models into a multiscale model using the DFBA approximation (Figure 1).

2.1. Construction and Modeling of the Quorum-Sensing Gene Regulatory Network

The following elements were used to construct and model the QS gene regulatory network, (i) bibliomic information about the genes involved in the QS process, (ii) a language/platform for the biological representation of the retrieved bibliomic data, and (iii) a language/platform for the mathematical modeling of the biological representation.
Bibliomic information was retrieved from generic and specialized databases, such as the Scopus and Web of Science bibliographic indexes and specific biological databases. In addition, several search queries were used to find relevant documents and database entries that provided the maximum number of chemical species for modeling the biological process.
The QS network model was built using the standard System Biology Markup Language (SBML) format. The CellDesigner 4.4 software [24] was used to model the biological processes by applying the Systems Biology Graphical Notation (SBGN) and to create the mathematical representation of the biological model with a MathML layer. The resulting model was based on bibliomic analysis and information stored in the GenomeNet (https://www.genome.jp, accessed on 5 March 2016), KEGG [25,26], BioCyc [27,28], UniProtKB/Swiss-Prot [29] and PubChem [30,31] databases, as well as on the search of annotated genes and their corresponding mRNAs and proteins to manually curate the model. Based on the biological model, a system of 114 ordinary differential equations (ODEs) representing the interactions between chemical species was established in a manually curated QS network model (Figure 2).
For QS modeling, the initial concentrations of genes, QS regulatory proteins, and available cytosolic Fe3+ were fixed (data available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1). In addition, each reaction incorporated the kinetic constant (k) values (a total of 103 kinetic parameters data are available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1) and the links of each interaction between related molecules. The corresponding ODEs for mRNAs and proteins in the transcription and translation processes are exemplified in Equations (1) and (2), respectively [32].
d Y m R N A dt =   k a   Y g e n e   k b   Y m R N A k c   Y m R N A
dYProtein dt =   k 1   Y m R N A   k 2   YProtein k 3 YProtein
where Y represents the chemical species, k a     the transcription rate constant, k b   and k 1 the translation rate constants, k c   and k 3   the degradation rate constants, and k 2   the consumption rate constant of each reaction in the system.
To construct the deterministic QS model, elementary kinetics, based on the law of mass action for all chemical species, was assumed for all reactions in the system. QS model simulations were run using SOSlib [33] in the CellDesigner software, which solves the rigidity problem of the ODEs and initiates numerical integration using the backward differentiation formula (BDF) or the Adams-Moulton (AM) method to calculate x(t) for a series of time points [34].

Simulation Scenario Conditions for the Quorum-Sensing Network

After testing the stability of the QS network by the behavior of signal molecules, intracellular ferribactin production, and intracellular and extracellular PVD production, a single condition was set for the different simulation scenarios, namely, the initial concentration of the extracellular signal molecule PQS (E-PQS). This ranged from 0.01 μM to 0.1 μM with 0.01-unit intervals and from 0.01 μM to 0.1 μM with 0.1-unit intervals for a total of 20 simulation scenarios (Table 1). These simulation scenarios were performed to emulate the behavior of in vitro cultures of P. aeruginosa. When the microbial population density increases in these cultures, the concentration of QS signal molecules that regulate the expression of bacterial phenotypes, such as PVD, also increases [9,35,36]. In addition, the simulations gave results for changing metabolite concentrations and reaction fluxes in time ( μ mol   s 1 L 1 ) . The latter ones were used in the subsequent modeling of the P. aeruginosa metabolic network.

2.2. Construction of the Pseudomonas aeruginosa Metabolic Network

The in silico construction of the P. aeruginosa metabolic network required a previously published generic cell model of the bacterium, the specific metabolic pathway of PVD synthesis, and software for modeling the hundreds of reaction equations representing the metabolic pathways.
The computational modeling was based on three data sources: (i) The genome-scale metabolic network model iMO1056 [37], which was to be used as a template, because it was the first genome-scale metabolic network construction of P. aeruginosa published in the literature; (ii) the P. aeruginosa genome annotation in the PseudoCAP database; and (iii) the metabolic pathway of PVD synthesis in the MetaCyc database (PseudoCyc), which was extended by bibliomics. As in the IMO1056 model, each reaction was manually mapped against data in biological databases (KEGG, BioCyc [27], ModelSEED [38], MetanetX [39], TransportDB [40], TCDB [41]) for preliminary curation. This process updated and complemented the network regarding stoichiometry, directionality, biological evidence, and subcellular location of each reaction, as well as in the gene-protein-reaction (GPR) associations. Furthermore, specific reactions involved in the PVD metabolic synthesis, some exchange reactions involving QS signal molecules, reactions for PVD transport, and reactions in bacterial culture (Luria-Bertani, LB) medium, were considered for computational modeling (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).

2.2.1. Curation of the Pseudomonas aeruginosa Metabolic Network

Curing a metabolic process refers to making the network mathematically feasible by removing pathologies such as metabolites that are not produced or consumed by any reaction in the model, eliminating thermodynamically infeasible cycles, and validating the production of metabolites in the biomass reaction.
Before curating the P. aeruginosa metabolic network, the list of reactions in the network was represented as a stoichiometric matrix ( S i j ) containing the ratio of the stoichiometric coefficient of each metabolite ( i ) in the reaction ( j ) . Next, this network was curated according to the following steps: (i) validation of the production of the biomass precursors by reverse engineering [42], (ii) addition of the PVD metabolite to the biomass reaction by increasing the value of the PVD coefficient until the metabolic pathway of PVD synthesis was activated, (iii) detection and solution of network pathologies such as root no-production and root no-consumption metabolites [43,44], and (iv) search and solution of thermodynamically infeasible cycles (TICs) in the network [45]. In this work, the proposal to add PVD to the biomass reaction was inspired by the work of Amara et al. [46], Prigent et al. [47], and Kim et al. [48], in which efforts have been made to optimize the production of secondary metabolites, such as virulence factors, in genome-scale metabolic models.

2.2.2. Modeling of the Pseudomonas aeruginosa Metabolic Network Using a Steady-State FBA Approximation

The P. aeruginosa metabolic network modeling using a steady-state FBA approximation [17] comprised the representation of a system of differential equations coupled to an objective function, i.e., the biomass maximization (an optimization problem, a term borrowed from engineering). The optimization problem was formulated as described in Equation (3) and solved as a linear programming (LP) problem using the CPLEX solver in the GAMS software (General Algebraic Modeling System: https://www.gams.com/, accessed on 20 March 2023).
M a x i m i z e     μ = v j B i o m a s s  
subject to:
j = 1 J   S i j   v j = 0 ,             i = 1 , , I
L B j   v j U B j ,               j = 1 , , J   j R   j I R J
v O 2 = K 1
v G l c = K 2
v A T P = K 3
where μ is the objective function, I   and J   represent the total number of metabolites and reactions, respectively;   S i j is the stoichiometric matrix for each metabolite i in the reaction j ; and v j is the flux of the reaction j expressed in mmol   gDW 1 h 1 ; and   L B j and U B j are the minimum and maximum fluxes that can adopt every reaction j   given by the thermodynamic feasibility of reactions. The fluxes for oxygen ( v O 2 ) and glucose ( v G l c ) uptake, and ATP ( v A T P ) production were limited at the rates of 10   mmol   gDW 1 h 1   ( K 1 and K 2 ), and 10   mmol   gDW 1 h 1   ( K 3 ) [37].

2.3. Combining the QS Gene Regulatory Network and the Pseudomonas aeruginosa Metabolic Network Models into an Integrative Model

The QS gene regulatory network and the P. aeruginosa metabolic network represent two types of biological processes; both are part of a living cell but differ in time scale, spatial location, control mechanisms, and chemical species that usually work separately. However, both networks share chemical species, such as the enzymes encoded by the QS network genes that catalyze some biochemical reactions in the P. aeruginosa metabolic network. The challenge was to combine these two networks to function as an integrated multiscale model. The work of Mallmann and collaborators inspired the methodological proposal to solve the challenge [49]. Thus, the challenge was solved by (i) Finding the shared reactions between the two models: the QS gene regulatory network model and Pseudomonas aeruginosa metabolic network model, and (ii) using the fluxes of the reactions shared by both model as constraints in the multi-stage FBA (see Section 2.3.2) and DFBA (see Section 2.3.3) simulations. The two networks were found to share nine reactions involved in the synthesis and transport of QS signal molecules (PQS, 3O-C12-HSL, and C4-HSL), ferribactin and PVD (Table 2).

2.3.1. Design of Simulation Scenarios

The integrated model was evaluated through several simulated scenarios expecting its results to match the known cellular behavior of an in vitro bacterial culture of P. aeruginosa [50,51]. For the simulations, six scenarios were selected for subsequent simulations using the multi-stage FBA and DFBA approximations. Scenarios Sc1, Sc2, and Sc3 included reactions involved in the synthesis of QS signal molecules and those involved in ferribactin synthesis (Sc1), PVD synthesis (Sc2), or PVD transport (Sc3). Scenarios Sc4, Sc5, and Sc6 included reactions involved in the transport of QS signal molecules and those involved in ferribactin synthesis (Sc4), PVD synthesis (Sc5), or PVD transport (Sc6) (Table 3). These simulation scenarios were performed to evaluate which combination of reactions in the multi-stage FBA (see Section 2.3.2) and DFBA (see Section 2.3.3) simulations showed a more significant change in the fluxes’ distribution. Scenario 6 (Sc6) was chosen.

2.3.2. Simulation Using the Multi-Stage FBA Approximation

Under the multi-stage FBA approximation, both network models were combined by fixing in the metabolic model, as equality constraints, the fluxes of the shared reactions, obtained from the QS network simulations according to each simulation scenario where E-PQS concentration was changed (Table 1). The units of the fluxes of the QS reactions obtained from the QS network simulations ( μ mol   s 1 L 1 ) were adjusted to mmol   gDW 1 h 1 [52]. Simulations were run at a time scale of 2000 intervals in the six simulation scenarios (Table 3). The flux value obtained from the QS model in each time interval of the shared reactions, according to the reactions involved in the different scenarios (Table 3), was fixed respectively for each time interval in the multi-stage FBA and DFBA (see Section 2.3.3) simulations. In addition, the optimization problem (Equation (3)) was modified for the multi-stage FBA (Equation (4)) and solved as a linear programming (LP) problem using the CPLEX v.12.6.0.0 solver in GAMS (https://www.gams.com/latest/docs/S_CPLEX.html, accessed on 20 March 2023).
M a x i m i z e     μ = v j B i o m a s s  
Subject to:
j = 1 J   S i j *   v j = 0 ,               i = 1 , , I              
L B j   v j U B ,                 j = 1 , , J   j R   j I R J                        
v O 2 = K 1          
v G l c = K 2          
v A T P   = K 3          
v j = v j ^               j ^ Q S   R e a c t i o n s
  t t 0 ,   t f
where j ^   is the subset of J -type reactions denoting the reactions shared by the QS and P. aeruginosa metabolic networks, w is the instance to evaluate μ in j ^ -type reactions, α is the final instance to evaluate μ in j ^ -type reactions and β is the last j ^ -type reaction to be fixed. For this model, the fluxes of the shared reactions ( v j ^ ) were fixed according to each simulation scenario (Table 3). Finally, scenario 6 (Sc6) showed a more significant change in the fluxes’ distribution.

2.3.3. Simulation Using the DFBA Approximation

The behavior of the integrated model over time was observed by dynamic simulation using the DFBA approximation. Thus, observing the influence of the QS network on the P. aeruginosa metabolic network and the interaction between both networks over time was possible. Furthermore, the DFBA approximation proposed by Mahadevan et al. [22] was used to solve the optimization problem shown in Equation (5).
M a x i m i z e j = 1 N c j · v j t
Subject to:
z i t + Δ t = z i t j = 1 J S i j · v j t · X t · Δ t             i   ϵ   1 , ,   I e x t r a c e l l u l a r
X t + Δ t = X t + µ · X t · Δ t
j = 1 J S i j · v j = 0             i   ϵ   1 , ,   I e x t r a c e l l u l a r
v j m i n < v j t < v j m a x               j   ϵ   1 , ,   J
c ^ z i t ,   v j t 0                                                                        
z i t 0               z i t 0 = z i , 0               i   ϵ   1 ,   ,   I e x t r a c e l l u l a r
X t 0           X t 0 = X 0                                                                    
Δ t = t f t 0 G               G   ϵ   0 G
v j = v j ^             j ^   ϵ   Q s   R e a c t i o n s
  t   ϵ   t 0 ,   t f      
where z i is the extracellular metabolite concentration, z i , 0 and X 0   are the initial conditions for each metabolite and biomass concentration, respectively, µ is the specific cell growth rate, c j is the reaction weight, c ^ z , v is the vector of nonlinear constraints (substrate consumption kinetics), t 0 and t f correspond to the initial and final times (0 and 19 h, respectively, at 0.0095-h intervals), G is the number of intervals used to discretize the time and j ^   is the subset of J -type reactions denoting the reactions shared by the QS and P. aeruginosa metabolic networks. Simulations were run to obtain (i) the flux profile of each reaction, (ii) the growth rate profile, (iii) the biomass concentration profile, and (iv) the concentration profiles over time of the metabolites of interest, i.e., the QS signal molecules, glucose, and PVD. The optimization problem was solved as a nonlinear programming problem (NLP) using the CONOPT v3.15N solver in GAMS software (https://www.gams.com/latest/docs/S_CONOPT.html, accessed on 20 March 2023). This model includes the equality and inequality constraints of the system (thermodynamics) and the fluxes of the QS network reactions fixed in the P. aeruginosa metabolic network (according to scenarios Table 3) as equality constraints according to the results of each simulation scenario where E-PQS concentration was changed (Table 1). The following conditions were set for the simulation: 0.0095-time step size and three orthogonal placement points for a Legendre polynomial; also, initial source concentrations of carbon (glucose = 55.5 mM), nitrogen (NH4+ = 104.127 mM), phosphate (Pi = 3.88 mM), sulfate (SO42− = 0.286 mM), biomass (X = 0.1 g/L) and oxygen (O2 = 1 mM). The initial concentrations of glucose, nitrogen, phosphate, and sulfate were calculated according to the approximate composition of the Luria Broth (LB) culture medium reported in the technical specifications of the commercial culture medium.

2.4. Cultures of Pseudomonas aeruginosa Strain PAO1

The integrated model, named CCBM1146, was evaluated qualitatively by comparing the behavior pattern of PVD production and biomass production between data obtained from in silico simulations and those obtained from in vitro cultures of P. aeruginosa strain PAO1 (collection of the research group Comunicación y comunidades bacterianas, Department of Biology, Faculty of Sciences, Universidad Nacional de Colombia, Bogotá) was grown in Luria Broth medium (10 g/L tryptone, 5 g/L yeast extract, 10 g/L NaCl) in a 1-L BioStat®A bioreactor (Sartorius Stedim Biotech, Goettingen, Germany), equipped with a disc-type impeller, a ring-type diffused bubble aeration system, pH and pO2 probes, a chiller for temperature control, and an aeration system for continuous and automatic control of the laboratory airflow. Cultures were grown in triplicate batches under controlled pH = 7, Tº = 37 °C, agitation = 200 rpm, oxygen saturation = 20%, and air flow = 2 L*min−1 (2 vvm). Cultures were incubated for 24 h, and aliquots were collected at 1-h intervals.

2.4.1. Evaluation of Bacterial Growth

Aliquots of 1 mL of the P. aeruginosa cultures were taken to measure the absorbance at 600 nm by spectrophotometry. The mean values were analyzed by logistic regression using the StatSoft statistical software to fit the data according to Equation (6) and obtain the following specific growth equation for strain PAO1.
v 2 = E x p   ( ( a + b ) * v 1 1 + E x p ( a + ( b ) * v 1 ) )
where v 2 corresponds to absorbance and v 1   to time. Values calculated for a ( 1.950216 )   and b   ( 0.234159 ) , showed an adequate adjustment of the data with R = 0.9915 and explained   variance = 98.316 % .

2.4.2. Evaluation of Biomass Production

After incubation was complete, three aliquots (1 mL) of the in vitro cultures of P. aeruginosa were taken and centrifuged at 14,000 rpm for 5 m at 4 °C. The pellets were discarded, and the supernatants were stored at −80 °C for subsequent determination of the PVD profile following the protocol of Meyer and Abdallah [53,54,55]. Total biomass was measured as dry weight. Briefly, at the end of the culture incubation, three aliquots (50 mL) of the cultures were collected and centrifuged at 14,000 rpm for 5 m. The pellets were washed with sterile water, centrifuged again at 14,000 rpm for 5 m, dried at 50 °C, transferred to a drying chamber, and weighed.

3. Results

3.1. The Quorum-Sensing Gene Regulatory Network for Pyoverdine Expression in Pseudomonas aeruginosa: A Deterministic Model

The QS network model described here consists of 114 subcellular chemical species, including proteins, small molecules, genes and their corresponding mRNAs, and nine complexes containing different molecules (data available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1). In addition, this QS network model comprised 103 biochemical reactions, including DNA transcription, mRNA translation, protein complex formation, inhibition reactions, and molecule-to-molecule interactions with positive, negative, or unknown effects (Figure 2). All reactions and interactions between them had the kinetic rate constant (k) values reported in the literature. The values of k not reported in the literature were assumed to be within the range of published k values for similar reactions (data available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1). For model construction, a system of 114 ODE was generated using a total of 103 kinetic parameters and represented using the SBGN graphical notation of the SBML models.
Figure 2. The quorum-sensing network regulates the expression of pyoverdine in Pseudomonas aeruginosa. The QS network was constructed by analysis of bibliomic information, schematized in CellDesigner 4.4, and modeled with the SBML ODE Solver Library (SOSlib). The chemical species are represented in SBML format: Genes by straight yellow squares; mRNAs by dark green diagonal squares; proteins by light green and brown oval squares; single molecules by bright green and pink ovals; and O2, Fe2+, and Fe3+ ions by blue circles. Gray and white squares represent complexes of simple molecules and macromolecular complexes, respectively. Red circles represent the degradation processes of all proteins included in the model. Arrow colors represent the following interactions between chemical species: Purple for transcription and translation processes; black for complex formation processes; yellow for complex cleavage; pink and green for the diffusion of signal molecules; red and blue for negative and positive regulatory processes, respectively; and orange for protein interactions in specific reactions. (Each symbol and its meaning are available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1).
Figure 2. The quorum-sensing network regulates the expression of pyoverdine in Pseudomonas aeruginosa. The QS network was constructed by analysis of bibliomic information, schematized in CellDesigner 4.4, and modeled with the SBML ODE Solver Library (SOSlib). The chemical species are represented in SBML format: Genes by straight yellow squares; mRNAs by dark green diagonal squares; proteins by light green and brown oval squares; single molecules by bright green and pink ovals; and O2, Fe2+, and Fe3+ ions by blue circles. Gray and white squares represent complexes of simple molecules and macromolecular complexes, respectively. Red circles represent the degradation processes of all proteins included in the model. Arrow colors represent the following interactions between chemical species: Purple for transcription and translation processes; black for complex formation processes; yellow for complex cleavage; pink and green for the diffusion of signal molecules; red and blue for negative and positive regulatory processes, respectively; and orange for protein interactions in specific reactions. (Each symbol and its meaning are available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1).
Metabolites 13 00659 g002
The QS network was modeled in the twenty simulation scenarios described in Table 1. The system was set up in all of these scenarios with an initial concentration of 1.0 μM of all QS genes and super-regulator proteins. In addition, extracellular Fe3+ availability and PVD transport system proteins were fixed to a starting concentration of 2.0 μM (data available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1). In all simulation scenarios, the concentration of the E-PQS signal molecule was modified to emulate the behavior of in vitro cultures of P. aeruginosa, where population growth increases as a function of time (more cells per unit volume) (Table 1). The performance of the QS network model was understood as synthesis of the chemical species of interest, i.e., intracellular production of the QS signal molecules (homoserine-lactones 3O-C12-HSL and C4-HSL, and quinolone PQS), cytosolic production of ferribactin, and periplasmic production of PVD (Figure 3).
The first simulation was run under the “initial conditions” scenario (Table 1). This simulation tested the intracellular production of the chemical species of interest (data available in Mendeley Data, https://doi.org/10.17632/2xzzkmnpfx.1). In this simulation, increasing concentrations of QS signal molecules were not considered because the purpose was to assess the basal capacity of the model to reproduce the QS circuits and produce QS signal molecules intracellularly. In silico intracellular production of the chemical species of interest corresponded to simulated concentrations of 0.0234 μM for 3O-C12-HSL, 0.0102 μM for C4-HSL, 0.00099 μM for PQS; 0.0158 μM for ferribactin, and 0.0635 μM for PVD (Figure 3). These values are relevant because they test whether the model behaves in agreement with data reported in the literature [56,57,58].
The low concentration of these chemical species showed that the system on its own produced a basal amount of signal molecules, ferribactin, and PVD in response to the concentration of Fe3+ available in the extracellular medium (simulated at an initial concentration of 2.0 μM). Subsequently, simulations run according to the scenarios in Table 1 showed that the higher the concentration of E-PQS, the higher the intracellular PVD production. The experimentation performed allowed the evaluation of the sensitivity of the output variable [intracellular PVD] to marginal changes in the input parameter [E-PQS] (Appendix A).
The highest PVD yield (0.122 μM) was achieved with 0.6 μM E-PQS (Sc16; ID = PQSE06 in Table 1). On the other hand, the highest concentrations of E-PQS (Sc17 to Sc20) were related to a decrease in PVD production (Figure 4A), indicating that the saturation point of the system was reached at 0.6 μM E-PQS. At the saturation point in the cell, in the particular case of the signaling molecule PQS, when there is a PQS “maximum” concentration inside the cell—which depends directly on the QS phenomenon [7,8]—no new extracellular PQS molecules can enter the cell, thus decreasing the virulence factors such as QS-regulated PVD production over time [13,14,15,16].
A dynamic equilibrium was reached in each simulation scenario, meaning that all chemical species in the model reached a constant value after some time interval. Contrary to the effect of the initial E-PQS concentration on intracellular PVD production, the final extracellular PVD concentration was not affected in any of the simulated scenarios (Figure 4B). However, the system evidenced a change in the trajectory of extracellular PVD production.

3.2. Pseudomonas aeruginosa Metabolic Network Model CCBM1146: An Improved Version of the Genome-Scale Metabolic Model iMO1056

The iMO1056 model, including its biomass reaction [37], was the starting point. First, the data associated with all reactions in the model were reviewed, followed by a reverse-engineering curation process. Next, 90 genes, 120 metabolites, and 41 reactions, including those in the metabolic pathway of PVD synthesis, were selected to fill the gaps in the model. Finally, PVD was added to the biomass reaction as a metabolite to obtain the improved model CCBM1146 (Table 4).
The curation process was complemented by (i) evaluation of metabolite production and consumption in the biomass reaction by reverse engineering; (ii) detection and resolution of pathologies in the 880 metabolites of the network: 29.51% had pathologies, 17.25% were resolved and 12.26% with unresolved pathologies remained in the model; (iii) detection and resolution of TICs: 23 TICs were identified and entirely resolved in the network. Finally, the proposed model, CCBM1146, comprised 1123 reactions, 880 metabolites, 136 protein complexes and transport proteins, and 1146 genes coding for 826 enzymes. Each reaction was manually assigned a metabolic system and its corresponding labeled subsystems according to the ontology pathway (https://biocyc.org/searchhelp.shtml#ontology_search, accessed on 15 October 2018) established by the BioCyc database (Figure 5). Analysis of the BioCyc metabolic system labels revealed that the CCBM1146 model consisted mainly of reactions of central metabolism, such as biosynthesis of cofactors, prosthetic groups, and electron transporters, followed by transport and exchange reactions involved in amino acid biosynthesis and degradation, lipopolysaccharide biosynthesis, and reactions involved in the biosynthesis of lipid molecules for cell wall formation (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).
The definitive CCBM1146 model and the steady-state FBA approximation were used to obtain the flux distribution of reactions in the P. aeruginosa metabolic network. The flux value (the rate at which a metabolite is formed as a function of the available biomass [mmol/gWD*h−1]) obtained for the objective function of the model was 0.55 h−1 (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1). After simulation, 352 reactions with active fluxes were obtained, including those for synthesis, transport, and exchange of the QS signal molecules (PQS, 3O-C12-HSL, and C4-HSL) and PVD. This relatively small number of reactions is principally due to the structure of the biomass equation. Also, it could be related to the addition of PVD to the biomass reaction, which triggers the activation of a small set of essential reactions of the central metabolism (including biomass metabolites and PVD synthesis reactions). Furthermore, since PVD is a secondary metabolite, these reactions result in PVD production upon activating necessary intermediary metabolism and secondary metabolism reactions in the model.

3.3. Integrative Model Simulations Evidenced the Influence of Quorum-Sensing Signaling on Pyoverdine Biosynthesis in Pseudomonas aeruginosa Cultures

As described in materials and methods, the design of the integrative model required the combination of the QS regulatory gene network model and the P. aeruginosa metabolic network model. Since these networks have different types of biological language, the flux values of nine reactions shared by both networks were considered to be useful for merging the behavior of chemical species in both types of network models. These nine shared reactions were involved in the synthesis and transport of QS signal molecules, ferribactin synthesis, and PVD synthesis and transport. After applying the multi-stage FBA approximation in six simulation scenarios of the integrated model (Table 3), the Sc3 scenario was chosen for subsequent simulations. This scenario evidenced the highest sensitivity of the CCBM1146 model to changes in the fluxes of reactions involved in the metabolic synthesis of the three QS signal molecules and transport of PVD. Subsequently, twenty simulations were run under the multi-stage FBA approximation in scenarios with varying initial E-PQS concentration (Table 1) and fixing the flux values of four of the reactions shared by the QS and the P. aeruginosa metabolic network models included in the Sc3 scenario (Table 3). Figure 6A shows a dynamic response of the integrative model objective function in the first four hours of in silico simulations; this response was related to changes in signal intensity derived from the QS model (deterministic model). This response was evident in all simulation conditions of the CCBM1146 model. Therefore, the sensitivity of the CCBM1146 model can be biologically interpreted as follows: the higher the E-PQS concentration at the beginning of bacterial growth, the higher the bacterial metabolic demand, and the greater the deceleration of the objective function (biomass maximization) (Figure 6A).
Twenty simulations run using the DFBA approximation and under conditions identical to those of the multi-stage FBA approximation showed similar results (Figure 6B), supporting the consistency of the proposed CCBM1146 model. During the first hours of simulation, the objective function showed a dynamic response related to changes in the signal intensity derived from the QS network in all simulation scenarios, i.e., as the QS signal increased, the behavior of the objective function changed. Subsequently, the objective function reached the same value in all simulation scenarios at approximately the 14-h interval. Finally, the value of the objective function decreased to zero simultaneously with the overall consumption of the available glucose (carbon source) in the model.
In addition, biomass profiles and concentrations of metabolites of interest were obtained under the DFBA approximation. Both in silico (Figure 7A) and in vitro (Figure 7B) data showed an exponential increase in biomass profiles up to hour 14 of culture of P. aeruginosa. A similar exponential increasing trend in biomass concentration was also observed over time up to the stationary phase. Therefore, it is possible to suggest that the biomass profile predicted (Figure 7A) by the CCBM1146 model properly represents the in vitro behavior of P. aeruginosa (Figure 7B).
In silico and in vitro PVD profiles were also compared (Figure 8). In silico, the PVD value increased in the first 8 h of simulated culture and reached a stable concentration around 9 h (Figure 8A). The in vitro PVD profile showed an increasing trend in PVD concentration over time. The maximum PVD concentration was reached around 10 h of culture. It coincided roughly with the mid-exponential phase of bacterial growth in vitro (Figure 8B) when signal molecules were produced. Bacteria responded to those signal molecules in the stationary phase [59].
However, a variation in PVD concentration was observed in vitro, probably related to typical PVD oxidation over time in a natural system. Moreover, the oxidation changed the characteristic green color of PVD in vitro due to the quinoline-based cyclic fluorescent chromophore, which is responsible for the bright fluorescence of pyoverdine [60].
In silico, glucose concentration varied over time in all simulation conditions (Figure 9A); the initial glucose concentration (55 mM) was consumed entirely at approximately 14 h. In addition, in silico profiles of the QS signal molecules, i.e., 3O-C12-HSL (Figure 9B), C4-HSL (Figure 9C), and PQS (Figure 9D), also showed their concentrations increased up to hour 10. This time roughly corresponded to the mid-exponential phase of P. aeruginosa growth in vitro and agreed with data from the literature, according to which cells have a basal production of signal molecules that reach a maximum concentration in the mid-exponential phase of growth.

4. Discussion

A thorough understanding of complex biological systems’ behavior requires studying and comprehending the interplay of several processes that often occur at different spatial and temporal scales in the biological system as a whole [1,2,61]. The biological problem addressed in this work was the synthesis of the virulence factor PVD in P. aeruginosa. The synthesis of PVD requires bacterial communication through the QS phenomenon to produce signal molecules that regulate the expression of pvd genes encoding enzymes that catalyze the synthesis of the virulence factor PVD. Three types of biological processes are involved in PVD synthesis, namely: (i) the QS signaling pathway, (ii) a regulatory pathway for pvd gene expression, and (iii) the effect of these two processes on the P. aeruginosa metabolic network for PVD biosynthesis. These three processes were encompassed here in a multiscale model to explain the influence of the QS phenomenon on the metabolic pathway of PVD biosynthesis in P. aeruginosa.

4.1. The QS Gene Regulatory Network Model Emulates the Natural Behavior of Pseudomonas aeruginosa

The deterministic in silico model of the QS gene regulatory network in P. aeruginosa proposed in this work is the first approximation that includes the three QS systems reported for this microorganism. It also incorporates the autoregulatory mechanisms of these three QS systems as well as the mechanisms responsible for the expression of pvd genes, the transport of PVD to the extracellular space, and the chelation and entry of iron into the bacterium through the PVD/Fe3+ complex. The data obtained from QS model simulations are supported by the fact that any system in dynamic equilibrium tends to reach a steady state, which translates into an equilibrium that resists external forces of change (perturbations). When the system is perturbed, the regulatory systems—in this case, the QS systems—respond to the signals emitted to establish a new equilibrium; this function is executed by the QS systems themselves and is known as feedback control. All biological processes that integrate and coordinate the functioning of living organisms are examples of the homeostatic regulation required by the system to ensure cell survival. Thus, the biological system is driven to the homeostatic state expected from its behavior, which is likely to maintain stability while adapting to the optimal conditions for survival. This argument is supported by the work of Gonzalez et al. [62], which states that the steady state is the inherent state of biological systems in the environment. Therefore, analyses of bacterial systems during the steady state have been the basis for most conclusions on bacterial network modeling [62].
The proposed deterministic QS network model emulated in silico the natural behavior of P. aeruginosa in culture, i.e., as the bacterial population increased (more cells per unit volume), the extracellular concentration of QS signal molecules per unit volume also increased. It was also evident that the signal molecules of all three QS systems were synthesized intracellularly and that as their concentration decreased, PVD production increased; that is, the QS network model reproduced known biological behavior and worked according to the expected cellular dynamics concerning the QS systems (Figure 3). It is important to note that a complex modeling of the internal regulation of the QS network in P. aeruginosa was achieved, as well as the possibility to simulate extracellular conditions that can lead to an increase in bacterial population density by modifying the concentration of the exogenous QS signal molecules. As shown in Figure 4A, the simulated intracellular PVD production was directly proportional to the extracellular PQS concentration (increase in QS strength) under the simulation conditions posed for each scenario (Table 1). However, after reaching a peak concentration, PVD production decreased, indicating that the system reached a saturation point with the following range of dynamic equilibrium. This result could be related to an increase in the intracellular Fe2+ concentration resulting from the autoregulatory mechanisms of the system. In the cytosol, the ferrous ion (Fe2+) forms a complex with the Fur protein, Fur/Fe2+, that binds to the promoter regions of pvd genes (iron-repressible or iron-regulated genes), thus inhibiting their transcription. Under limiting intracellular iron concentrations, the Fur-mediated repression decreases, and a positive transcriptional regulation of pvd genes ensues [63,64,65]. The proposed deterministic QS network model reproduces this regulatory phenomenon that is directly related to the production of PVD, the formation of the PVD/Fe3+ complex at the extracellular level, its subsequent entry into the cell via an ABC-type transport system, and finally, to the release of Fe2+ into the cytosolic space.
On the other hand, it should be noted that the simulated extracellular PVD (Figure 4B) behaved differently from intracellular PVD in response to the increase of extracellular PQS levels. In the first few hours, the system showed a change in the trajectory of the extracellular PVD concentration. This different dynamic response could be attributed to the PVD transport mechanism involving the surface receptor FpvA, the anti-sigma factor FpvR [66,67], and the protein responsible for energy transduction to import molecules through the TonB protein located in the outer membrane of the bacterium [7,8,14,15,67,68,69]. According to Arevalo-Ferro et al., the expression of TonB is also regulated by QS in P. aeruginosa [70]. Thus, the results of simulations performed with the QS Network can be supported, as mentioned previously, on the basis that any system in dynamic equilibrium tends to reach a stable state, which translates into an equilibrium that resists external forces of change.
In addition, the model acquired the emerging property of resilience, characteristic of adaptive biological systems. Although this property was not modeled, it became apparent when the model was disturbed as the QS regulatory system responded to outputs to establish a new dynamic equilibrium that reached a steady state. The generated equilibrium resisted disturbances, understood as external forces of change. The in silico results were similar to those reported on the natural behavior of QS mechanisms in P. aeruginosa [67,71].

4.2. The Proposed Integrative CCBM1146 Model Helps to Infer the Influence of the QS Phenomenon on the PVD Metabolic Biosynthesis in Pseudomonas aeruginosa

This study proposed a systemic approach to design an integrative multiscale model to understand the influence of the QS phenomenon on the expression of the PVD metabolic phenotype of P. aeruginosa by combining two networks with different temporal and spatial scales, namely, the QS communication network responsible for regulating the expression of pvd genes and the bacterial metabolic network of PVD biosynthesis.
The CCBM1146 metabolic network model is an improved version of the genome-scale metabolic network iMO1056 developed by Oberhardt et al. [37]. The CCBM1146 model includes reactions of the central metabolism, all reactions for the biosynthesis of PVD and QS signal molecules, as well as those involved in their corresponding transport and exchange. The strategy used to combine the two networks was based on obtaining from the QS network simulations the fluxes of the reactions shared by the two networks to set them in the corresponding metabolic reactions in the CCBM1146 model as constraints of the optimization problem in the DFBA. The network combination strategy was evaluated by simulations in the integrative CCBM1146 model under a multi-stage optimization approach to obtain the values of changes in the objective function (biomass maximization) over time without considering the changes in metabolite concentrations. Subsequently, the model was run under a DFBA approximation to obtain the distribution of reaction fluxes over time and the concentration profiles of biomass and metabolites of interest. As shown in Figure 6A, the results of the multi-stage FBA modeling evidenced a different dynamic response of the objective function related to changes in the signal strength from the QS network model in all simulation scenarios. As the QS signal increased, there was a noticeable change in behavior of the objective function, which decreased to a minimum and then increased to a stable value in all simulation scenarios. The decrease in the objective function value (Figure 6A) could be attributed to its degeneracy by one of the elements involved in the optimization, which does not support its value to increase. Degeneracy also occurs when the constraint in the optimization problem increases [72,73]. According to the study by Wintermute et al., “the FBA generally cannot predict a unique rate for all fluxes. A solution that maximizes growth rate is typically mathematically degenerate, describing a region in flux space rather than a single point. Solution degeneracy is a well-described problem in systems like metabolism, which are flexible, internally redundant, and underdetermined by data” [72]. Alternatively, it may also occur because the availability of resources does not imply that the model immediately exploits them in favor of maximizing the objective function. Instead, it can be interpreted as a period of adjustment or adaptation of all variables to optimize the maximization of the objective function. However, as the influence of the constraint increases—understood as a more significant influence of the QS signal—the point at which the value of the objective function does not vary significantly over time (steady state) is reached in a shorter time.
The DFBA results were similar to those from the multi-stage FBA approximation (Figure 6B). However, under the DFBA approximation, the objective function reached a stable value around the hour 14 and then dropped to zero in all simulation scenarios. The point at which the objective function reached a zero value coincided with the point of depletion of the carbon source (glucose). This fact could explain the results, considering that once the total glucose consumption is reached, the maximization of the objective function cannot take a value other than zero. In in vitro cultures, bacterial growth reaches the stationary phase when the nutrients in the medium—especially the carbon source—have been consumed, indicating that the bacteria can no longer grow. Growth curves of Pseudomonas aeruginosa in cultures incubated under controlled conditions—intended to validate some of the results of the CCBM1146 model—showed that the stationary phase was reached between 14 h and 15 h of culture (Figure 10), similar to the in silico point at which the value of the objective function dropped to zero.
The divergences observed between the in vitro and in silico results may be due to several factors, including the growth conditions of each system. For the in silico model, glucose was considered the sole carbon source; its concentration, as well as those for nitrogen, phosphate, and sulfate sources, were estimated from the composition of the LB culture medium. The LB culture medium is a complex mixture of nutrients because its main constituents are yeast extract and tryptone, so the concentration of each nutrient could not be accurately determined. In contrast, in the in vitro model, the bacteria grew in cultures prepared with LB medium. There, they could take advantage of the total concentration of nutrients and utilize carbon sources other than glucose, which may be reflected in bacterial growth. Furthermore, it is important to note that the in vitro growth conditions were optimal for the microorganism, even better than those in vivo.
Finally, it should be noted that in the in vitro cultures, the bacterial population produced the QS signal by itself; i.e., the culture medium was not supplemented with synthetic QS signal molecules. Nevertheless, the pattern of behavior of the biomass profile (Figure 7B), and PVD concentration profile (Figure 8B) were adequately reproduced by the in silico results obtained in the integrative model under the DFBA approximation (Figure 7A and Figure 8A).
The methodological strategy employed in this study and the results indicate that the proposed multiscale CCBM1146 model offers valuable insights into the influence of the QS phenomenon on the metabolism and subsequent behavior of P. aeruginosa and other microorganisms related to the synthesis of virulence factors such as PVD. In addition, the CCBM1146 model could help to infer the metabolic behavior of P. aeruginosa that arise in the presence of different concentrations of QS signal molecules. Therefore, the present study also can be a starting point for comprehending the biosynthesis of other siderophores regulated by QS in P. aeruginosa, such as Pyochelin (PCH), which, like PVD, has a high iron affinity, and the newly discovered narrow-spectrum metallophore Pseudopaline (PSP), that acts as a virulent factor too and is involved in nickel and zinc uptake [6,75,76]. This study made possible the design of an original methodology based on tools from different areas of knowledge, such as biological network reconstruction methods from systems biology, genomics, and bibliomics, as well as engineering methods, such as optimization, from the area of logistics and operation research; together, they provided a way to study the hierarchies of biological systems in a unified way, according to the principles of systems biology.

5. Conclusions

From a systemic perspective, this work presents a methodological strategy to design a novel multi-scale and multi-class model proposal by combining two classes of models with different scales. The CCBM1146 model helped to characterize and explain the complex and emerging behavior derived from the interactions between these two models, which would have been impossible by studying each model or scale separately. Moreover, the CCBM1146 model served to understand and infer the natural influence of PQS intensity on the metabolic PVD biosynthesis in P. aeruginosa. In the future, this integrative model could help to infer the different metabolic phenotypes that may arise among other microorganisms when exposed to different concentrations of the signal molecules of their own QS circuits responsible for regulating the synthesis of their virulence factors. Furthermore, the DFBA approximation applied to the CCBM1146 model evidenced its descriptive capability for the profiles of biomass and PVD, showing a similar behavioral pattern of an in vitro culture of P. aeruginosa; this makes it more valuable and suitable for evaluating, in the future and from holistic and dynamic perspectives, how different conditions can affect the behavior of P. aeruginosa for the synthesis of different virulence factors regulated by QS.
Moreover, the model could be used to test experimental conditions in vitro. Finally, this work proposes, for the first time, the integration of the quorum-sensing gene regulatory network with the P. aeruginosa metabolic network for PVD biosynthesis. Combining these two network models was possible by fixing the fluxes of reactions shared by both models as system constraints in the multi-stage FBA and DFBA approximation. Thus, it was possible to model the influence of the QS phenomenon on the P. aeruginosa metabolism to PVD biosynthesis as a function of QS signal intensity.
However, though in this work a sensitivity analysis was performed for the QS network model, for further modeling works under the steady-state and dynamic-state flux balance analysis approach of the CCBM1146 integrated model, it is recommended that a sensitivity analysis be performed.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/metabo13050659/s1. File S1: Data Descriptor Information.

Author Contributions

Conceptualization, D.C.C.-B., C.A.-F. and A.F.G.B.; Methodology; D.C.C.-B. and A.F.G.B.; Software: D.C.C.-B. and A.F.G.B.; Validation, D.C.C.-B.; Formal Analysis, D.C.C.-B.; Investigation, D.C.C.-B.; Resources, D.C.C.-B. and C.A.-F. (P. aeruginosa strain PAO1); Data Curation, D.C.C.-B.; Writing—Original Draft Preparation, D.C.C.-B.; Writing—Review & Editing, D.C.C.-B. and A.F.G.B.; Visualization, D.C.C.-B.; Supervision, C.A.-F. and A.F.G.B.; Project Administration, D.C.C.-B. and A.F.G.B.; Funding Acquisition, D.C.C.-B. All authors have read and agreed to the published version of the manuscript.

Funding

This project was funded by the Colombian Ministry of Science and Technology (Colciencias) through Grant ID: 747 of 2015 and by the OMICAS program: In silico Multiscale Optimization of Sustainable Agricultural Crops (Infrastructure and Validation in Rice and Sugarcane) (Optimización Multiescala In silico de Cultivos Agrícolas Sostenibles) (Infraestructura y validación en Arroz y Caña de Azúcar), sponsored within the Colombian Scientific Ecosystem by The World Bank, COLCIENCIAS, ICETEX, the Colombian Ministry of Education and the Colombian Ministry of Industry and Tourism under Grant ID: FP44842-217-2018 and OMICAS Award ID: 792-61187. This funding institutions had no role in the study design, data collection, analysis, manuscript preparation, or publication decision.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Quorum-Sensing Model for the Pyoverdine Expression in P. aeruginosa. https://doi.org/10.17632/2xzzkmnpfx.1 (https://data.mendeley.com/datasets/2xzzkmnpfx), and P. aeruginosa Genome-scale Metabolic Network-CCBM1146. https://doi.org/10.17632/y9htx3fcjm.1 (https://data.mendeley.com/datasets/y9htx3fcjm). Details data can be available in the Supplementary Material.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

2C7-3OH-4(1H)Q: 2-heptyl-3-hydroxy-4-quinolone. 3O-C12-HSL: N-(3-Oxododecanoyl)-L-homoserine lactone. AM: Adams-Moulton method. BDF: Backward differentiation formula. C4-HSL: N-Butanoyl-L-homoserine lactone. CCBM1146: Genome-scale metabolic model of Pseudomonas aeruginosa proposed in this work. DFBA: Dynamic flux balance analysis. FBA: Flux balance analysis. GAMS: General Algebraic Modeling System software. GPR: Gene-Protein-Reaction. LB: Luria-Bertani culture medium. LP: Linear programming problem. NLP: Non-linear programming problem. NRPS: Non-ribosomal peptide synthase. ODEs: Ordinary differential equations system. PAO1: Pseudomonas aeruginosa strain PAO1. E-PQS: Extracellular Pseudomona quinolone signal. PQS: Pseudomona quinolone signal. PVD: Pyoverdine. PVD-Fe3+: Ferri-pyoverdine complex. QS: Quorum-sensing. SBGN: Systems Biology Graphical Notation. SBML: Systems Biology Markup Language. TICs: Thermodynamically infeasible cycles.

Appendix A

Sensitivity analysis of the QS network model. This analysis evaluated the sensitivity of the output variable [intracellular PVD] to marginal changes in the input parameter [E-PQS].
Table A1. Sensitivity analysis of the QS network model by subsets.
Table A1. Sensitivity analysis of the QS network model by subsets.
SubsetInput Parameter Value [E-PQS]Maximum Input Parameter Value in Subset [E-PQS]Output Variable Value [PVD]Average Output Variable Value in Subset [PVD]
Sensitivity Subset 1 = 0.68110.1201779860.118943557
0.90.120450180
0.80.120766334
0.70.121045165
0.60.121118254
0.50.120687054
0.40.119372306
0.30.116454998
0.20.110419737
Sensitivity Subset 2 = 0.520.10.10.0971730800.082671936
0.090.095075669
0.080.092768928
0.070.090219930
0.060.087397210
0.050.084268237
0.040.080803399
0.030.076979188
0.020.072803560
0.010.068338620
The value in the input parameter box corresponds to the value of the E-PQS concentration used in each experiment in Table 1. The value in the output variable box corresponds to the value of the maximum intracellular PVD concentration obtained in each experiment of Table 1 and is plotted in Figure 4A. The bold value in the input parameter box indicates the value of the parameter at which the output variable acquires its maximum value in each subset. [PVD] indicates the intracellular pyoverdine concentration and is the output variable. [E-PQS] indicates the extracellular Pseudomona quinolone signal and is the input parameter.
The sensitivity of each subset box was calculated according to the following equation:
S e n s i t i v i t y S u b s e t = m a x .   p a r m i n .   p a r 2 m a x .   o b s .   p a r   A v e r a g e   v a r i a b l e m a x . v a r i a b l e
where m a x .   p a r is the maximum input parameter value in the subset, m i n .   p a r is the minimum input parameter value in the subset, m a x .   o b s .   p a r is the input parameter value with which the maximum output variable value is observed, A v e r a g e   v a r i a b l e   is the average value of the output variable in the subset, and m a x .   v a r i a b l e   is the maximum output variable value in the subset. Thus, the numerator of the equation measures the parameter’s proportion in the subset, and the denominator of the equation represents the variable’s proportion in the subset. The sensitivity value obtained for subset 1 was 0.68, which means that, on average, each element in subset 1 of the parameters represents 68% of the variable. The sensitivity value obtained for subset 2 was 0.52, which means that, on average, each element of subset 2 of parameters represents 52% of the variable.
The difference in sensitivity values between subset 1 (0.68) and subset 2 (0.52) indicates that the [PVD] output variable is more sensitive to changes in [E-PQS] input parameter in the range of 1.0 to 0.1 than in the range of 0.09 to 0.01. This suggests that the system may be more susceptible to changes in the [E-PQS] input parameter in subset 1. In addition, each element in subset 1 represents 68% of the output variable, which means that, on average, the system may be more susceptible to changes in [E-PQS] input parameter in the range of 1.0 to 0.1, suggesting that the system might be more sensitive to changes in the [E-PQS] input parameter in subset 1. It also means that, on average, a change in an [E-PQS] input parameter value in subset 1 (range 1.0 to 0.1) has a more significant impact on the [PVD] output variable compared to a change in an [E-PQS] input parameter value in subset 2 (range 0.09 to 0.01). In other words, the system is more sensitive to changes in [E-PQS] input parameter in the range of values in subset 1 than in the range of subset 2. Sensitivity analysis on different subsets of data can provide information about the robustness of the model to variations in the input parameters. For example, if the sensitivity is high in a specific range of [E-PQS] input parameter values, this could indicate that the model is less robust to changes in [E-PQS] input parameter in that range.
Table A2. Rate of change to each input parameter value [E-PQS].
Table A2. Rate of change to each input parameter value [E-PQS].
Input Parameter Value [E-PQS]Output Variable Value [PVD]Rate of Change to Each Input Parameter Value [E-PQS]
10.1201779868.32
0.90.1204501807.47
0.80.1207663346.62
0.70.1210451655.78
0.60.1211182544.95
0.50.1206870544.14
0.40.1193723063.35
0.30.1164549982.58
0.20.1104197371.81
0.10.0971730801.03
0.090.0950756690.95
0.080.0927689280.86
0.070.0902199300.78
0.060.0873972100.69
0.050.0842682370.59
0.040.0808033990.50
0.030.0769791880.39
0.020.0728035600.27
0.010.0683386200.15
The value in the input parameter box corresponds to the value of the E-PQS concentration used in each experiment in Table 1. The value in the output variable box corresponds to the value of the maximum intracellular PVD concentration obtained in each experiment in Table 1 and is plotted in Figure 4A. [PVD] indicates the intracellular pyoverdine concentration and is the output variable. [E-PQS] indicates the extracellular Pseudomona Quinolone Signal and is the input parameter.
The rate of change of each input parameter box was calculated according to the following equation:
R a t e   o f   C h a n g e = I p v   O v v
where I p v is the adopted [E-PQS] input parameter value in each experiment, and O v v is the resultant value of the [PVD] output variable in each experiment; this [PVD] output variable value depends on the [E-PQS] input parameter value assigned for each experiment. This rate-of-change analysis provides information on how the output variable [PVD] responds to changes in the input parameter [E-PQS] over different ranges of values. As the value of the input parameter [E-PQS] decreases, the rate of change also decreases. This suggests that the system response (PVD) to changes in the input parameter [E-PQS] is more pronounced for higher values of the input parameter [E-PQS]. The relationship between the input parameter [E-PQS] and the rate of change is not linear. In the range of input parameter [E-PQS] from 1.0 to 0.2, the rate of change decreases faster at input parameter [E-PQS] from 0.2 to 0.01. This indicates that the relationship between input parameter [E-PQS] and output variable [PVD] may have nonlinear behavior at different ranges of input parameter [E-PQS].
In summary, this rate-of-change analysis provides additional information on how the output variable [PVD] responds to changes in the input parameter [E-PQS] at different ranges of values.

References

  1. Martins, M.L.; Ferreira, S.C.; Vilela, M.J. Multiscale models for biological systems. Curr. Opin. Colloid Interface Sci. 2010, 15, 18–23. [Google Scholar] [CrossRef]
  2. Walpole, J.; Papin, J.A.; Peirce, S.M. Multiscale computational models of complex biological systems. Annu. Rev. Biomed. Eng. 2013, 15, 137–154. [Google Scholar] [CrossRef] [PubMed]
  3. Warner, H.V.; Sivakumar, N.; Peirce, S.M.; Lazzara, M.J. Multiscale computational models of cancer. Curr. Opin. Biomed. Eng. 2019, 11, 137–144. [Google Scholar] [CrossRef]
  4. Thiele, I.; Fleming, R.M.T.; Que, R.; Bordbar, A.; Diep, D.; Palsson, B.O. Multiscale Modeling of Metabolism and Macromolecular Synthesis in E. coli and Its Application to the Evolution of Codon Usage. PLoS ONE 2012, 7, e45635. [Google Scholar] [CrossRef]
  5. Biggs, M.B.; Papin, J.A. Novel Multiscale Modeling Tool Applied to Pseudomonas aeruginosa Biofilm Formation. PLoS ONE 2013, 8, e78011. [Google Scholar] [CrossRef]
  6. Ghssein, G.; Ezzeddine, Z. A Review of Pseudomonas aeruginosa Metallophores: Pyoverdine, Pyochelin and Pseudopaline. Biology 2022, 11, 1711. [Google Scholar] [CrossRef]
  7. Visca, P.; Imperi, F.; Lamont, I.L. Pyoverdine siderophores: From biogenesis to biosignificance. Trends Microbiol. 2007, 15, 22–30. [Google Scholar] [CrossRef]
  8. Visca, P.; Imperi, F.; Lamont, I.L. Pyoverdine Synthesis and its Regulation in Fluorescent Pseudomonads. In Microbial Siderophores; Springer: Berlin/Heidelberg, Germany, 2007; Volume 12, pp. 135–163. [Google Scholar]
  9. Givskov, M.; Rasmussen, T.B.; Ren, D.; Balaban, N. Bacterial Cell-to-cell Communication (Quorum Sensing). In Control of Biofilm Infections by Signal Manipulation; Springer: Berlin/Heidelberg, Germany, 2008; pp. 13–38. [Google Scholar]
  10. Lazdunski, A.M.; Ventre, I.; Bleves, S. Cell–Cell Communication: Quorum Sensing and Regulatory Circuits in Pseudomonas aeruginosa. In Pseudomonas; Ramos, J.-L., Filloux, A., Eds.; Springer: Dordrecht, The Netherlands, 2007; pp. 279–310. [Google Scholar]
  11. Imperi, F.; Visca, P. Subcellular localization of the pyoverdine biogenesis machinery of Pseudomonas aeruginosa: A membrane-associated “siderosome”. FEBS Lett. 2013, 587, 3387–3391. [Google Scholar] [CrossRef]
  12. Ringel, M.T.; Dräger, G.; Brüser, T. PvdN Enzyme Catalyzes a Periplasmic Pyoverdine Modification. J. Biol. Chem. 2016, 291, 23929–23938. [Google Scholar] [CrossRef]
  13. Cobessi, D.; Celia, H.; Folschweiller, N.; Schalk, I.J.; Abdallah, M.A.; Pattus, F. The crystal structure of the pyoverdine outer membrane receptor FpvA from Pseudomonas aeruginosa at 3.6 Å resolution. J. Mol. Biol. 2005, 347, 121–134. [Google Scholar] [CrossRef]
  14. Greenwald, J.; Hoegy, F.; Nader, M.; Journet, L.; Mislin, G.L.; Graumann, P.L.; Schalk, I.J. Real time fluorescent resonance energy transfer visualization of ferric pyoverdine uptake in Pseudomonas aeruginosa. A role for ferrous iron. J. Biol. Chem. 2007, 282, 2987–2995. [Google Scholar] [CrossRef] [PubMed]
  15. Voulhoux, R.; Filloux, A.; Schalk, I.J. Pyoverdine-mediated iron uptake in Pseudomonas aeruginosa: The Tat system is required for PvdN but not for FpvA transport. J. Bacteriol. 2006, 188, 3317–3323. [Google Scholar] [CrossRef] [PubMed]
  16. Beare, P.A.; For, R.J.; Martin, L.W.; Lamont, I.L. Siderophore-mediated cell signalling in Pseudomonas aeruginosa: Divergent pathways regulate virulence factor production and siderophore receptor synthesis. Mol. Microbiol. 2003, 47, 195–207. [Google Scholar] [CrossRef] [PubMed]
  17. Orth, J.D.; Thiele, I.; Palsson, B.Ø. What is flux balance analysis? Nat. Biotechnol. 2010, 28, 245–248. [Google Scholar] [CrossRef] [PubMed]
  18. Raman, K.; Chandra, N. Flux balance analysis of biological systems: Applications and challenges. Brief. Bioinform. 2009, 10, 435–449. [Google Scholar] [CrossRef]
  19. Gianchandani, E.P.; Chavali, A.K.; Papin, J.A. The application of flux balance analysis in systems biology. Wiley Interdiscip. Rev. Syst. Biol. Med. 2010, 2, 372–382. [Google Scholar] [CrossRef]
  20. Maranas, C.D.; Zomorrodi, A.R. Optimization Methods in Metabolic Networks; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2016. [Google Scholar]
  21. Serrano-Bermúdez, L.M.; Barrios, A.F.G.; Maranas, C.D.; Montoya, D. Clostridium butyricum maximizes growth while minimizing enzyme usage and ATP production: Metabolic flux distribution of a strain cultured in glycerol. BMC Syst. Biol. 2017, 11, 58. [Google Scholar] [CrossRef]
  22. Mahadevan, R.; Edwards, J.S.; Doyle, F. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys. J. 2002, 83, 1331–1340. [Google Scholar] [CrossRef]
  23. Serrano-Bermúdez, L.M.; Barrios, A.F.G.; Montoya, D. Clostridium butyricum population balance model: Predicting dynamic metabolic flux distributions using an objective function related to extracellular glycerol content. PLoS ONE 2018, 13, e0209447. [Google Scholar] [CrossRef]
  24. Funahashi, A.; Matsuoka, Y.; Jouraku, A.; Kitano, H.; Kikuchi, N. Celldesigner: A Modeling Tool for Biochemical Networks. In Proceedings of the 2006 Winter Simulation Conference, Monterey, CA, USA, 3–6 December 2006; pp. 1707–1712. [Google Scholar]
  25. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef]
  26. Kanehisa, M.; Furumichi, M.; Tanabe, M.; Sato, Y.; Morishima, K. KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 2017, 45, D353–D361. [Google Scholar] [CrossRef]
  27. Karp, P.D.; Billington, R.; Caspi, R.; Fulcher, C.A.; Latendresse, M.; Kothari, A.; Keseler, I.M.; Krummenacker, M.; Midford, P.E.; Ong, Q.; et al. The BioCyc collection of microbial genomes and metabolic pathways. Brief. Bioinform. 2019, 20, 1085–1093. [Google Scholar] [CrossRef]
  28. Caspi, R.; Altman, T.; Billington, R.; Dreher, K.; Foerster, H.; Fulcher, C.A.; Holland, T.A.; Keseler, I.M.; Kothari, A.; Kubo, A.; et al. The MetaCyc database of metabolic pathways and enzymes and the BioCyc collection of Pathway/Genome Databases. Nucleic Acids Res. 2014, 42, 473–479. [Google Scholar] [CrossRef] [PubMed]
  29. Bairoch, A.; Apweiler, R. The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res. 2000, 28, 45. [Google Scholar] [CrossRef] [PubMed]
  30. Kim, S.; Thiessen, P.A.; Bolton, E.E.; Chen, J.; Fu, G.; Gindulyte, A.; Han, L.; He, J.; He, S.; Shoemaker, B.A.; et al. PubChem Substance and Compound databases. Nucleic Acids Res. 2016, 44, D1202–D1213. [Google Scholar] [CrossRef] [PubMed]
  31. Kim, S.; Chen, J.; Cheng, T.; Gindulyte, A.; He, J.; He, S.; Li, Q.; Shoemaker, B.A.; Thiessen, P.A.; Yu, B.; et al. PubChem in 2021: New data content and improved web interfaces. Nucleic Acids Res. 2021, 49, D1388–D1395. [Google Scholar] [CrossRef] [PubMed]
  32. Aldridge, B.B.; Burke, J.M.; Lauffenburger, D.A.; Sorger, P.K. Physicochemical modelling of cell signalling pathways. Nat. Cell Biol. 2006, 8, 1195–1203. [Google Scholar] [CrossRef]
  33. Rainer, M.; Christoph, F.; Andrew, F.; Stefan, M.; James, L.; Akira, F. SoSLib; Github: Tokio, Japan, 2017. [Google Scholar]
  34. Machne, R.; Finney, A.; Muller, S.; Lu, J.; Widder, S.; Flamm, C. The SBML ODE Solver Library: A native API for symbolic and fast numerical analysis of reaction networks. Bioinformatics 2006, 22, 1406–1407. [Google Scholar] [CrossRef]
  35. Fuqua, W.C.; Winans, S.C.; Greenberg, E.P. Quorum sensing in bacteria: The LuxR-LuxI family of cell density-responsive transcriptional regulators. J. Bacteriol. 1994, 176, 269–275. [Google Scholar] [CrossRef]
  36. Fuqua, C.; Parsek, M.R.; Greenberg, E.P. Regulation of gene expression by cell-to-cell communication: Acyl-homoserine lactone quorum sensing. Annu. Rev. Genet. 2001, 35, 439–468. [Google Scholar] [CrossRef]
  37. Oberhardt, M.A.; Puchalka, J.; Fryer, K.E.; Santos, V.A.P.M.D.; Papin, J.A. Genome-Scale Metabolic Network Analysis of the Opportunistic Pathogen Pseudomonas aeruginosa PAO1. J. Bacteriol. 2008, 190, 2790–2803. [Google Scholar] [CrossRef] [PubMed]
  38. Seaver, S.M.D.; Liu, F.; Zhang, Q.; Jeffryes, J.; Faria, J.P.; Edirisinghe, J.N.; Mundy, M.; Chia, N.; Noor, E.; Beber, M.E.; et al. The ModelSEED Biochemistry Database for the integration of metabolic annotations and the reconstruction, comparison and analysis of metabolic models for plants, fungi and microbes. Nucleic Acids Res. 2021, 49, D575–D588. [Google Scholar] [CrossRef] [PubMed]
  39. Moretti, S.; Tran, V.D.T.; Mehl, F.; Ibberson, M.; Pagni, M. MetaNetX/MNXref: Unified namespace for metabolites and biochemical reactions in the context of metabolic models. Nucleic Acids Res. 2021, 49, D570–D574. [Google Scholar] [CrossRef] [PubMed]
  40. Elbourne, L.D.H.; Tetu, S.G.; Hassan, K.A.; Paulsen, I.T. TransportDB 2.0: A database for exploring membrane transporters in sequenced genomes from all domains of life. Nucleic Acids Res. 2017, 45, D320–D324. [Google Scholar] [CrossRef]
  41. Saier, M.H., Jr.; Reddy, V.S.; Moreno-Hagelsieb, G.; Hendargo, K.J.; Zhang, Y.; Iddamsetty, V.; Lam, K.J.K.; Tian, N.; Russum, S.; Wang, J.; et al. The Transporter Classification Database (TCDB): 2021 update. Nucleic Acids Res. 2021, 49, D461–D467. [Google Scholar] [CrossRef]
  42. Senger, R.S.; Papoutsakis, E.T. Genome-scale model for Clostridium acetobutylicum: Part I. Metabolic network resolution and analysis. Biotechnol. Bioeng. 2008, 101, 1036–1052. [Google Scholar] [CrossRef]
  43. Kumar, V.S.; Dasika, M.S.; Maranas, C.D. Optimization based automated curation of metabolic reconstructions. BMC Bioinform. 2007, 8, 212. [Google Scholar]
  44. Saadat, N.P.; van Aalst, M.; Ebenhöh, O. Network Reconstruction and Modelling Made Reproducible with moped. Metabolites 2022, 12, 275. [Google Scholar] [CrossRef]
  45. Schellenberger, J.; Lewis, N.E.; Palsson, B.Ø. Elimination of Thermodynamically Infeasible Loops in Steady-State Metabolic Models. Biophys. J. 2011, 100, 544–553. [Google Scholar] [CrossRef]
  46. Amara, A.; Takano, E.; Breitling, R. Development and validation of an updated computational model of Streptomyces coelicolor primary and secondary metabolism. BMC Genom. 2018, 19, 519. [Google Scholar] [CrossRef]
  47. Prigent, S.; Nielsen, J.C.; Frisvad, J.C.; Nielsen, J. Reconstruction of 24 Penicillium genome-scale metabolic models shows diversity based on their secondary metabolism. Biotechnol. Bioeng. 2018, 115, 2604–2612. [Google Scholar] [CrossRef] [PubMed]
  48. Kim, H.U.; Charusanti, P.; Lee, S.Y.; Weber, T. Metabolic engineering with systems biology tools to optimize production of prokaryotic secondary metabolites. Nat. Prod. Rep. 2016, 33, 933–941. [Google Scholar] [CrossRef]
  49. Mallmann, J.; Heckmann, D.; Bräutigam, A.; Lercher, M.J.; Weber, A.P.; Westhoff, P.; Gowik, U. The role of photorespiration during the evolution of C4 photosynthesis in the genus Flaveria. Elife 2014, 3, e02478. [Google Scholar] [CrossRef] [PubMed]
  50. Monod, J. The Growth of Bacterial Cultures. Annu. Rev. Microbiol. 1949, 3, 371–394. [Google Scholar] [CrossRef]
  51. Folsom, J.P.; Richards, L.; Pitts, B.; Roe, F.; Ehrlich, G.D.; Parker, A.; Mazurie, A.; Stewart, P.S. Physiology of Pseudomonas aeruginosa in biofilms as revealed by transcriptome analysis. BMC Microbiol. 2010, 10, 294. [Google Scholar] [CrossRef]
  52. Milo, R.; Phillips, R. Cell Biology by the Numbers, 1st ed.; Garland Science: New York, NY, USA, 2015. [Google Scholar]
  53. Fallahzadeh, V.; Ahmadzadeh, M.; Sharifi, R. Growth and pyoverdine production kinetics of Pseudomonas aeruginosa 7NSK2 in an experimental fermentor. J. Agric. Technol. 2010, 6, 107–115. [Google Scholar]
  54. Imperi, F.; Tiburzi, F.; Visca, P. Molecular basis of pyoverdine siderophore recycling in Pseudomonas aeruginosa. Proc. Natl. Acad. Sci. USA 2009, 106, 20440–20445. [Google Scholar] [CrossRef]
  55. Meyer, J.M.; Abdallah, M.A. The Fluorescent Pigment of Pseudomonas fluorescens: Biosynthesis, Purification and Physicochemical Properties. J. Gen. Microbiol. 1978, 107, 319–328. [Google Scholar] [CrossRef]
  56. Pearson, J.P.; Gray, K.M.; Passador, L.; Tucker, K.D.; Eberhard, A.; Iglewski, B.H.; Greenberg, E.P. Structure of the autoinducer required for expression of Pseudomonas aeruginosa virulence genes [density-dependent transcription/gene activation/las genes/N-acyl homoserine lactone/3-oxo-N-(tetrahydro-2-oxo-3-furanyl)dodecanamide. J. Biochem. 1994, 91, 197–201. [Google Scholar]
  57. Venturi, V. Regulation of quorum sensing in Pseudomonas. FEMS Microbiol. Rev. 2006, 30, 274–291. [Google Scholar] [CrossRef]
  58. Venturi, V.; Weisbeek, P.; Koster, M. Gene regulation of siderophore-mediated iron acquisition in Pseudomonas: Not only the Fur repressor. Mol. Microbiol. 1995, 17, 603–610. [Google Scholar] [CrossRef] [PubMed]
  59. Gilbert, D.; Heiner, M.; Ghanbar, L.; Chodak, J. Spatial quorum sensing modelling using coloured hybrid Petri nets and simulative model checking. BMC Bioinform. 2019, 20, 173. [Google Scholar] [CrossRef] [PubMed]
  60. Nadal-Jimenez, P.; Koch, G.; Reis, C.R.; Muntendam, R.; Raj, H.; Jeronimus-Stratingh, C.M.; Cool, R.H.; Quax, W.J. PvdP is a tyrosinase that drives maturation of the pyoverdine chromophore in Pseudomonas aeruginosa. J. Bacteriol. 2014, 196, 2681–2690. [Google Scholar] [CrossRef]
  61. Twycross, J.; Band, L.R.; Bennett, M.J.; King, J.R.; Krasnogor, N. Stochastic and deterministic multiscale models for systems biology: An auxin-transport case study. BMC Syst. Biol. 2010, 4, 34. [Google Scholar] [CrossRef] [PubMed]
  62. Barrios, A.F.G.; Achenie, L.E.K. Escherichia coli autoinducer-2 uptake network does not display hysteretic behavior but AI-2 synthesis rate controls transient bifurcation. Biosystems 2010, 99, 17–26. [Google Scholar] [CrossRef]
  63. Leoni, L.; Orsi, N.; De Lorenzo, V.; Visca, P. Functional Analysis of PvdS, an Iron Starvation Sigma Factor of Pseudomonas aeruginosa Functional Analysis of PvdS, an Iron Starvation Sigma Factor of Pseudomonas aeruginosa. California Institu. J. Bacteriol. 2000, 182, 1481–1491. Available online: http://jb.asm.org/ (accessed on 17 September 2013). [CrossRef] [PubMed]
  64. Moon, C.D.; Zhang, X.-X.; Matthijs, S.; Schäfer, M.; Budzikiewicz, H.; Rainey, P.B. Genomic, genetic and structural analysis of pyoverdine-mediated iron acquisition in the plant growth-promoting bacterium Pseudomonas fluorescens SBW25. BMC Microbiol. 2008, 8, 7. [Google Scholar] [CrossRef]
  65. Vasil, M.L.; Ochsner, U.A.; Johnson, Z.; Colmer, J.A.; Hamood, A.N. The Fur-regulated gene encoding the alternative sigma factor PvdS is required for iron-dependent expression of the LysR-type regulator PtxR in pseudomonas aeruginosa. J. Bacteriol. 1998, 180, 6784–6788. [Google Scholar] [CrossRef]
  66. Diggle, S.P.; Matthijs, S.; Wright, V.J.; Fletcher, M.P.; Chhabra, S.R.; Lamont, I.L.; Kong, X.; Hider, R.C.; Cornelis, P.; Cámara, M.; et al. The Pseudomonas aeruginosa 4-Quinolone Signal Molecules HHQ and PQS Play Multifunctional Roles in Quorum Sensing and Iron Entrapment. Chem. Biol. 2007, 14, 87–96. [Google Scholar] [CrossRef]
  67. Jimenez, P.N.; Koch, G.; Thompson, J.A.; Xavier, K.B.; Cool, R.H.; Quax, W.J. The multiple signaling systems regulating virulence in Pseudomonas aeruginosa. Microbiol. Mol. Biol. Rev. 2012, 76, 46–65. [Google Scholar] [CrossRef]
  68. Chincholkar, S.B.; Chaudhari, B.L.; Rane, M.R. Microbial Siderophore: A State of Art. In Microbial Siderophores; Varma, A., Chincholkar, S.B., Eds.; Springer: Berlin/Heidelberg, Germany, 2007; pp. 233–242. [Google Scholar]
  69. Das, A.; Prasad, R.; Srivastava, A.; Giang, P.H.; Bhatnagar, K.; Varma, A. Fungal Siderophores: Structure, Functions and Regulation. In Microbial Siderophores; Springer: Berlin/Heidelberg, Germany, 2007; pp. 1–42. [Google Scholar]
  70. Arevalo-Ferro, C.; Hentzer, M.; Reil, G.; Görg, A.; Kjelleberg, S.; Givskov, M.; Riedel, K.; Eberl, L. Identification of quorum-sensing regulated proteins in the opportunistic pathogen Pseudomonas aeruginosa by proteomics. Environ. Microbiol. 2003, 5, 1350–1369. [Google Scholar] [CrossRef] [PubMed]
  71. Williams, P.; Cámara, M. Quorum sensing and environmental adaptation in Pseudomonas aeruginosa: A tale of regulatory networks and multifunctional signal molecules. Curr. Opin. Microbiol. 2009, 12, 182–191. [Google Scholar] [CrossRef] [PubMed]
  72. Wintermute, E.H.; Lieberman, T.D.; Silver, P.A. An objective function exploiting suboptimal solutions in metabolic networks. BMC Syst. Biol. 2013, 7, 98. [Google Scholar] [CrossRef] [PubMed]
  73. Schultz, A.; Qutub, A.A. Predicting internal cell fluxes at sub-optimal growth. BMC Syst. Biol. 2015, 9, 18. [Google Scholar] [CrossRef]
  74. Kim, D.; Chung, S.; Lee, S.; Choi, J. Relation of microbial biomass to counting units for Pseudomonas aeruginosa. Afr. J. Microbiol. Res. 2012, 6, 4620–4622. [Google Scholar]
  75. Ghssein, G.; Matar, S. Chelating Mechanisms of Transition Metals by Bacterial Metallophores ‘Pseudopaline and Staphylopine’: A Quantum Chemical Assessment. Computation 2018, 6, 56. [Google Scholar] [CrossRef]
  76. Lhospice, S.; Gomez, N.O.; Ouerdane, L.; Brutesco, C.; Ghssein, G.; Hajjar, C.; Liratni, A.; Wang, S.; Richaud, P.; Bleves, S.; et al. Pseudomonas aeruginosa zinc uptake in chelating environment is primarily mediated by the metallophore pseudopaline. Sci. Rep. 2017, 7, 17132. [Google Scholar] [CrossRef]
Figure 1. Methodological workflow for designing a computational model to infer the influence of the quorum-sensing phenomenon on the PVD metabolic synthesis in Pseudomonas aeruginosa from a systemic perspective. Stage 1. Construction of the QS gene regulatory network under a deterministic approach using a kinetics based on the mass action law. Stage 2. Construction of the P. aeruginosa metabolic network using the flux balance analysis approach. Stage 3. Integration of the QS network and the P. aeruginosa metabolic network into a multiscale model to be proposed and modeled under a dynamic flux balance analysis approach.
Figure 1. Methodological workflow for designing a computational model to infer the influence of the quorum-sensing phenomenon on the PVD metabolic synthesis in Pseudomonas aeruginosa from a systemic perspective. Stage 1. Construction of the QS gene regulatory network under a deterministic approach using a kinetics based on the mass action law. Stage 2. Construction of the P. aeruginosa metabolic network using the flux balance analysis approach. Stage 3. Integration of the QS network and the P. aeruginosa metabolic network into a multiscale model to be proposed and modeled under a dynamic flux balance analysis approach.
Metabolites 13 00659 g001
Figure 3. Simulation of the quorum sensing gene regulatory network model of Pseudomonas aeruginosa. Intracellular production of QS signal molecules (3O-C12-HSL, C4-HSL, PQS), ferribactin, and PVD under the “initial conditions” scenario (Table 1).
Figure 3. Simulation of the quorum sensing gene regulatory network model of Pseudomonas aeruginosa. Intracellular production of QS signal molecules (3O-C12-HSL, C4-HSL, PQS), ferribactin, and PVD under the “initial conditions” scenario (Table 1).
Metabolites 13 00659 g003
Figure 4. Simulated (A) intracellular and (B) extracellular pyoverdine production in Pseudomonas aeruginosa cultures varied as a function of initial E-PQS concentration. Simulations were run in 20 scenarios at varying initial concentrations of extracellular PQS (from 0.01 μM to 0.1 μM with 0.01-unit intervals and from 0.1 μM to 1.0 μM with 0.1-unit intervals (Table 1).
Figure 4. Simulated (A) intracellular and (B) extracellular pyoverdine production in Pseudomonas aeruginosa cultures varied as a function of initial E-PQS concentration. Simulations were run in 20 scenarios at varying initial concentrations of extracellular PQS (from 0.01 μM to 0.1 μM with 0.01-unit intervals and from 0.1 μM to 1.0 μM with 0.1-unit intervals (Table 1).
Metabolites 13 00659 g004aMetabolites 13 00659 g004b
Figure 5. Pseudomonas aeruginosa metabolic network model CCBM1146 for pyoverdine biosynthesis involves reactions of the central metabolism and secondary metabolism. The most representative metabolic systems corresponded to the biosynthesis of cofactors, prosthetic groups, and electron transporters, followed by reactions responsible for metabolite transport and exchange, and then amino acid biosynthesis and degradation. Another group of reactions was involved in the synthesis of lipopolysaccharides and other lipid molecules for cell wall formation. The reactions in the model were classified according to the metabolic system to which they are ascribed in the BioCyc database (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).
Figure 5. Pseudomonas aeruginosa metabolic network model CCBM1146 for pyoverdine biosynthesis involves reactions of the central metabolism and secondary metabolism. The most representative metabolic systems corresponded to the biosynthesis of cofactors, prosthetic groups, and electron transporters, followed by reactions responsible for metabolite transport and exchange, and then amino acid biosynthesis and degradation. Another group of reactions was involved in the synthesis of lipopolysaccharides and other lipid molecules for cell wall formation. The reactions in the model were classified according to the metabolic system to which they are ascribed in the BioCyc database (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).
Metabolites 13 00659 g005
Figure 6. In silico profile of the CCBM1146 model objective function under (A) the multi-stage FBA approximation and (B) the dynamic FBA approximation. Simulation results under the multi-stage FBA approximation were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixing the flux values of four reactions shared by the QS and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3).
Figure 6. In silico profile of the CCBM1146 model objective function under (A) the multi-stage FBA approximation and (B) the dynamic FBA approximation. Simulation results under the multi-stage FBA approximation were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixing the flux values of four reactions shared by the QS and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3).
Metabolites 13 00659 g006
Figure 7. In silico and In vitro biomass profiles of Pseudomonas aeruginosa cultures. (A). In silico biomass profile of P. aeruginosa cultures in the CCBM1146 model under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3). (B). In vitro biomass profile of P. aeruginosa cultures. This biomass profile was generated from average dry weight data of P. aeruginosa cultures grown under controlled laboratory conditions.
Figure 7. In silico and In vitro biomass profiles of Pseudomonas aeruginosa cultures. (A). In silico biomass profile of P. aeruginosa cultures in the CCBM1146 model under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3). (B). In vitro biomass profile of P. aeruginosa cultures. This biomass profile was generated from average dry weight data of P. aeruginosa cultures grown under controlled laboratory conditions.
Metabolites 13 00659 g007aMetabolites 13 00659 g007b
Figure 8. In silico and In vitro profiles of pyoverdine production in Pseudomonas aeruginosa cultures. (A). In silico profile of pyoverdine production in the CCBM1146 model under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3). (B). In vitro profile of pyoverdine production. This profile was generated from the average absorbance (450 nm) data [53,55] of cultures of P. aeruginosa grown under controlled laboratory conditions.
Figure 8. In silico and In vitro profiles of pyoverdine production in Pseudomonas aeruginosa cultures. (A). In silico profile of pyoverdine production in the CCBM1146 model under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3). (B). In vitro profile of pyoverdine production. This profile was generated from the average absorbance (450 nm) data [53,55] of cultures of P. aeruginosa grown under controlled laboratory conditions.
Metabolites 13 00659 g008
Figure 9. In silico (A) glucose profile, (B) QS signal molecule 3O-C12-HSL profile, (C) QS signal molecule C4-HSL profile, and (D) extracellular QS signal molecule PQS in Pseudomonas aeruginosa under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3).
Figure 9. In silico (A) glucose profile, (B) QS signal molecule 3O-C12-HSL profile, (C) QS signal molecule C4-HSL profile, and (D) extracellular QS signal molecule PQS in Pseudomonas aeruginosa under the DFBA approximation. Simulation results were generated in GAMS for the CCBM1146 model at varying initial E-PQS concentrations (Table 1) and fixed flux values of four reactions shared by the QS network and P. aeruginosa metabolic network models in the Sc3 scenario (Table 3).
Metabolites 13 00659 g009
Figure 10. In vitro growth profile of Pseudomonas aeruginosa. The growth curve shows an exponential phase starting at around hour 6 and a stationary phase between 19 h and 20 h, in agreement with Kim et al. [74]. The growth curve was generated from mean OD (600 nm) data of P. aeruginosa cultures (n = 3) that were analyzed by logistic regression. In the model equation, v2 represents the fitted absorbance (y-axis, dependent variable), and v1 (x-axis) represents time.
Figure 10. In vitro growth profile of Pseudomonas aeruginosa. The growth curve shows an exponential phase starting at around hour 6 and a stationary phase between 19 h and 20 h, in agreement with Kim et al. [74]. The growth curve was generated from mean OD (600 nm) data of P. aeruginosa cultures (n = 3) that were analyzed by logistic regression. In the model equation, v2 represents the fitted absorbance (y-axis, dependent variable), and v1 (x-axis) represents time.
Metabolites 13 00659 g010
Table 1. Simulation scenario conditions for the quorum-sensing network model.
Table 1. Simulation scenario conditions for the quorum-sensing network model.
Simulation ScenarioE-PQS [μM]ID Simulation ScenarioSimulation ScenarioE-PQS [μM]ID Simulation Scenario
Sc10.00Initial conditionsSc110.1PQSE01
Sc20.01PQSE001Sc120.2PQSE02
Sc30.02PQSE002Sc130.3PQSE03
Sc40.03PQSE003Sc140.4PQSE04
Sc50.04PQSE004Sc150.5PQSE05
Sc60.05PQSE005Sc160.6PQSE06
Sc70.06PQSE006Sc170.7PQSE07
Sc80.07PQSE007Sc180.8PQSE08
Sc90.08PQSE008Sc190.9PQSE09
Sc100.09PQSE009Sc201.0PQSE10
Table 2. Reactions shared by the Quorum-Sensing and Pseudomonas aeruginosa metabolic network models.
Table 2. Reactions shared by the Quorum-Sensing and Pseudomonas aeruginosa metabolic network models.
Quorum-Sensing
Network
Metabolic
Network
EC/TC NumberReaction Equation in the Metabolic Network Model
3O-C12-HSL production3O-C12-HSL synthesisEC 2.3.1.184[c]: 3oxddACP + amet <==> 5mta + ACP + h + n3oxdd-hsl
C4-HSL productionC4-HSL synthesisEC 2.3.1.184[c]: amet + butACP <==> 5mta + ACP + h + nb-hsl
PQS productionPQS
synthesis
EC 1.14.13.182[c]: fadh2 + h + hhq + o2 --> nad + h2o + pqs
C4-HSL
diffusion
C4-HSL transport-nb-hsl[c] <==> nb-hsl[e]
3O-C12-HSL
diffusion
3O-C12-HSL transport-n3oxdd-hsl[c] <==> n3oxdd-hsl[e]
PQS
diffusion
PQS transport-pqs[c] <==> pqs[e]
Ferribactin
production
Ferribactin synthesisEC 6.3.2.[c]: glu-L + tyr-L + (2) ser-L + arg-L + 24dab + (2) fohorn + lys-L + (2) thr-L --> fbn + (12) h2o + (2) h
PVD
production
PVD
synthesis
EC 1.14.18.[c]: fbn + o2 --> pvd1 + h2o
PVD
export
PVD
transport
TC-1.B.14.1.6pvd1[c] --> pvd1[e]
Table 3. Biochemical reactions involved in the integrative model simulation scenarios using FBA and DFBA approximations. Sc: Scenario number; PVD: pyoverdine.
Table 3. Biochemical reactions involved in the integrative model simulation scenarios using FBA and DFBA approximations. Sc: Scenario number; PVD: pyoverdine.
Metabolic ReactionSc1Sc2Sc3Sc4Sc5Sc6
3O-C12-HSL synthesisXXX
C4-HSL synthesisXXX
PQS synthesisXXX
C4-HSL transport XXX
3O-C12-HSL transport XXX
PQS transport XXX
Ferribactin synthesisX X
PVD synthesis X X
PVD transport X X
Table 4. Comparison of Pseudomonas aeruginosa metabolic network models iMO1056 and CCBM1146 (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).
Table 4. Comparison of Pseudomonas aeruginosa metabolic network models iMO1056 and CCBM1146 (data available in Mendeley Data, https://doi.org/10.17632/y9htx3fcjm.1).
Reactions and ComponentsModel
iMO1056CCBM1146
Metabolic reactions728774
Transport reactions150146
Biomass reaction11
Maintenance reaction11
Exchange reactions118120
Reactions for metabolite input from the culture medium8481
Total reactions10821123
Total metabolites760880
Total genes10561146
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Clavijo-Buriticá, D.C.; Arévalo-Ferro, C.; González Barrios, A.F. A Holistic Approach from Systems Biology Reveals the Direct Influence of the Quorum-Sensing Phenomenon on Pseudomonas aeruginosa Metabolism to Pyoverdine Biosynthesis. Metabolites 2023, 13, 659. https://doi.org/10.3390/metabo13050659

AMA Style

Clavijo-Buriticá DC, Arévalo-Ferro C, González Barrios AF. A Holistic Approach from Systems Biology Reveals the Direct Influence of the Quorum-Sensing Phenomenon on Pseudomonas aeruginosa Metabolism to Pyoverdine Biosynthesis. Metabolites. 2023; 13(5):659. https://doi.org/10.3390/metabo13050659

Chicago/Turabian Style

Clavijo-Buriticá, Diana Carolina, Catalina Arévalo-Ferro, and Andrés Fernando González Barrios. 2023. "A Holistic Approach from Systems Biology Reveals the Direct Influence of the Quorum-Sensing Phenomenon on Pseudomonas aeruginosa Metabolism to Pyoverdine Biosynthesis" Metabolites 13, no. 5: 659. https://doi.org/10.3390/metabo13050659

APA Style

Clavijo-Buriticá, D. C., Arévalo-Ferro, C., & González Barrios, A. F. (2023). A Holistic Approach from Systems Biology Reveals the Direct Influence of the Quorum-Sensing Phenomenon on Pseudomonas aeruginosa Metabolism to Pyoverdine Biosynthesis. Metabolites, 13(5), 659. https://doi.org/10.3390/metabo13050659

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop