Construction of spatial transcriptome-guided PA14 metabolic models

A key feature of MiMICS allows users to input multiple -omics-guided metabolic models, which may differ in predicted intracellular and extracellular metabolic parameter values, into a multi-scale metabolic framework (Fig 1). A set of four transcriptome-guided PA14 metabolic models was constructed by integrating an established PA14 GENRE [12] with a published spatial transcriptomic P. aeruginosa PA14 biofilm dataset [5] using the established algorithm RIPTiDe [11]. In addition, the PA14 GENRE was constrained on synthetic cystic fibrosis sputum medium (SCFM), replicating SCFM used in the experiment [5]. To avoid the computationally expensive generation of transcriptome-guided metabolic models for all ~292,000 PA14 biofilm cells measured, transcriptomic data representative of four unique metabolic states was extracted using the UMAP Leiden clustering method [5] and manual selection (refer to Methods). This approach generated four unique transcriptome-guided PA14 metabolic models. Each represented a unique metabolic state within the ten-hour PA14 biofilm, designated as (1) aerobic respiration, (2) denitrification, (3) denitrification + nitric oxide (NO) secretion, and (4) oxidative stress.

The aerobic respiration metabolic model, integrated with transcriptomic data of cells with high expression levels of ccoN1, which encodes the primary aerobic respiration oxidase [27], predicted high levels of oxygen uptake flux (Fig 2D). A denitrification PA14 metabolic model was generated from cells with high expression levels of genes encoding denitrification reductases, narG, napA, nirS, norB, and nosZ (Fig 2A), which convert nitrate to nitrogen to perform respiration in low oxygen environments [22] (complete PA14 GENRE denitrification pathway shown in Fig 2C). The denitrification metabolic model accurately predicted flux through the denitrification reductase reactions, resulting in characteristic nitrate uptake and nitrogen secretion fluxes (Fig 2B and 2D) [28]. A denitrification +NO secretion PA14 metabolic model was generated from cells with high narG and nirS expression, respectively encoding nitrate and nitrite reductase, but low norB expression, encoding NO reductase. This expression profile, considered to be limited by the norB expression, resulted in a metabolic model state with nitrate uptake and NO secretion flux (Fig 2D). When simulated by individual agents in MiMICS, this NO secretion rate was passed from the agent’s intracellular scale (GENRE) to extracellular scale (reaction-diffusion model), which is an essential step to produce the NO biofilm microenvironment that induces oxidative stress. Demonstrating the impact of transcriptomic data integration, the PA14 metabolic model without transcriptomic data integration, termed transcriptome-free in Fig 2D, predicted oxygen and nitrate uptake, but did not predict NO secretion denitrification metabolic processes (Fig 2D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. MiMICS incorporated four unique spatial-transcriptome-guided PA14 metabolic model states.

(A) Heatmap of experimental gene expression values of each PA14 metabolic model state. (B) Heatmap of predicted flux values of gene-encoded reactions for each PA14 metabolic model state. Flux values reported are absolute flux values normalized to the maximum absolute flux value among metabolic model states. (C) PA14 GENRE denitrification and oxidative stress pathway updated from iPau21 GENRE with NO secretion pathway, NO iron cytotoxic reaction, and katA-encoded oxidative stress reaction. (D) Predicted fluxes of metabolites exchanged with the extracellular environment for a transcriptome-free PA14 metabolic model and the four transcriptome-guided PA14 metabolic model states. Fluxes reported are normalized to the maximum flux among metabolic models.

https://doi.org/10.1371/journal.pcbi.1012031.g002

Lastly, an oxidative stress PA14 metabolic model was generated from cells with high expression of the gene katA, which encodes the antioxidant enzyme Catalase A [23] (Fig 2A). Accordingly, the PA14 GENRE was updated with NO cytotoxic and antioxidant Catalase A intracellular reactions [23]. Briefly, to mitigate cytotoxic NO binding to intracellular iron, Catalase A binds to NO, acting as a NO sink supplementary to NO degradation by NorB [23] (Fig 2C). Correspondingly, the oxidative stress metabolic model predicted NO was consumed (Fig 2D) and preferentially degraded by NorB for biomass synthesis before being diverted through the antioxidant Catalase A reaction to mitigate NO toxicity [23] (S2 Fig).

Constrained on replete SCFM, the aerobic PA14 metabolic model state predicted a biomass growth rate similar to an experimental growth rate of aqueous P. aeruginosa in SCFM [29] (S2 Fig). The denitrification +/-NO secretion PA14 metabolic models predicted the lowest biomass growth rates indicative of low oxygen conditions [30] (S2 Fig). In addition to the predicted NO secretion rates, the predicted oxygen, nitrate, and glucose exchange flux rates (Fig 2D) were provided as inputs to update the corresponding metabolite concentrations in the reaction-diffusion model. Due to the differences observed in gltB expression, which encodes a glucose binding protein, among metabolic model states (Figs 2A and Fig 2B), and previous reports relating denitrification metabolism and glucose uptake [31,32], glucose was chosen as the carbon source to simulate in the reaction-diffusion model.

Mechanistic incorporation of transcriptome-guided PA14 metabolic models into MiMICS

MiMICS employs ABM rules for an agent to select a metabolic model state from a set of transcriptome-guided metabolic models input by the user. In this study, a combination of mechanistic ABM rules informed from literature and the spatial transcriptomic dataset were implemented into MiMICS (Table 1). P. aeruginosa has been observed to upregulate genes related to aerobic and denitrification metabolism in aerobic and deplete oxygen conditions, respectively [22]. In addition, in the presence of extracellular NO, P. aeruginosa increases expression of katA, encoding the antioxidant Catalase A [23]. Thus, to decide between the four PA14 metabolic model states, agents compared their local oxygen and NO concentrations to respective concentration thresholds, [O₂]_t and [NO]_t (Table 1). The parameter value for [NO]_t was obtained from literature as the extracellular NO concentration that induced katA expression [23]. The parameter value for [O₂]_t was fit to experimental outputs (S4 Fig). An oxic O₂ threshold, 0.21 mM, had the smallest MiMICS model error, which suggests complete oxygen depletion was not essential to induce denitrification metabolism. Indeed, previous studies observed P. aeruginosa used denitrification in microaerobic conditions (~0.05 mM oxygen) as a possible supplementary or competitive respiration strategy to aerobic respiration [28]. The parameterized [O₂]_t value in the oxic range can be expected because the height (~10 μm) of the ten-hour biofilm was not large enough to predict significant oxygen depletion [33]. In the case that lower oxygen concentrations were present in the biofilm, other mechanisms that deplete oxygen near the biofilm, such as an oxygen boundary layer [34] or oxygen consumption by the planktonic phase above the biofilm may be present. These mechanisms serve as updates to simulate metabolite transport in upcoming MiMICS versions. Altogether, MiMICS suggests an oxygen transition point emerged within the PA14 biofilm that induced denitrification metabolism, but oxygen depletion mechanisms in regions outside of the biofilm may have been present.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Mechanistic ABM rules for P. aeruginosa agents to choose a PA14 metabolic model state.

P. aeruginosa agents compared their local oxygen and nitric oxide concentration to respective metabolite thresholds. A stochastic parameter, R_n, determined the agent’s decision between the denitrification +/- NO secretion metabolic models.

https://doi.org/10.1371/journal.pcbi.1012031.t001

A stochastic parameter, R_n, generated only for agents in the denitrification-inducing low oxygen and low NO condition, was used for agents to select between the two denitrification metabolic model states, one with and one without predicted NO secretion flux (Table 1). The R_n threshold was informed from experimental data, which suggested 6% of denitrification cells exhibited the norB-limiting gene expression profile associated with NO secretion (S2 Fig). This R_n parameter suggests cells stochastically expressed nirS or norB, which encode nitrite reductase and NO reductase respectively. Indeed, stochastic expression of denitrification genes has been observed in other bacterial species [35]. In addition, experimental observations of P. aeruginosa biofilms after ten hours showed high expression levels of pilA (S1 Fig), which encodes type IV pili protein (PilA) that facilitates surface motility and shapes biofilm structure [5,36]. As pili synthesis reactions were not in the current PA14 GENRE [37], simple surface motility ABM rules were incorporated to recapitulate the PA14 biofilm structure and total cell count (S3 Fig).

Transcriptome-guided MiMICS predicted microscale metabolic heterogeneity and NO microenvironment in PA14 biofilm

MiMICS was simulated with the mechanistic ABM rules controlling agent execution of a transcriptome-guided PA14 metabolic model state (Table 1) (referred to as transcriptome-guided MiMICS simulation). Qualitatively, in comparison to the experiment, the transcriptome-guided MiMICS simulation accurately predicted microscale, spatially confined biofilm niches in which cells heterogeneously upregulated denitrification and oxidative stress metabolic processes (Fig 3). In the simulation, these niches were located near the center of the biofilm colony, where microaerobic and variable NO concentrations were predicted (Fig 3). Thus, MiMICS simulations suggested that variable NO signal concentrations within microaerobic biofilm regions resulted in a heterogeneous population of cells using denitrification and oxidative stress metabolism co-existing within the same microscale niche. Qualitatively, experimental PA14 biofilms observed greater spatial dispersion of cells expressing denitrification and oxidative stress genes (Figs 3 and S5), suggesting other mechanisms may regulate denitrification gene expression, such as stochastic expression of the denitrification gene narG encoding nitrate reductase [35].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Transcriptome-guided MiMICS improved predictions of microscale metabolic heterogeneity and NO microenvironment in PA14 biofilm.

Shown are representative 3D renderings of PA14 biofilms grown for ten-hours from the experiment, and transcriptome-free and transcriptome-guided MiMICS simulations. Cells are colored according to their metabolic state. Shown are 2D yz slices (x = 160 μm, 22 μm, 100 μm for experimental, transcriptome-guided MiMICS, and transcriptome-free MiMICS, respectively) of cell metabolic states, and predicted oxygen and NO concentrations. The x-values were chosen to compare similar biofilm colony structures across experimental and simulation conditions. Experimental data was reconstructed from Dar and co-workers.

https://doi.org/10.1371/journal.pcbi.1012031.g003

To demonstrate the advantages of incorporating multiple transcriptome-guided metabolic models, MiMICS was simulated with the PA14 GENRE unconstrained by transcriptomic data (referred to as transcriptome-free MiMICS), which is the standard practice for current multi-scale metabolic frameworks [8,16,17]. Transcriptome-free MiMICS inaccurately predicted a homogeneous biofilm population with active flux through all intracellular denitrification metabolic reactions (Figs 3 and S5). In addition, transcriptome-free MiMICS did not predict extracellular NO in the biofilm microenvironment, resulting in agents lacking flux through the NO-induced oxidative stress reaction encoded by katA (Figs 3 and S5).

Transcriptome-guided MiMICS improved prediction accuracy of heterogenous PA14 biofilm metabolism

Quantitative validation of MiMICS was first performed by comparing the experimental percentage of total cells expressing a gene to the simulation’s percentage of total agents with active flux through metabolic reactions encoded by that respective gene (Fig 4). Compared to transcriptome-free MiMICS that predicted 100% of cells with active denitrification reaction flux, transcriptome-guided MiMICS accurately predicted reduced cell populations with active flux through denitrification narG-, napA-, nirS-, norB-, and nosZ-encoded reactions (Fig 4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Transcriptome-guided MiMICS improved prediction accuracy of heterogeneous intracellular metabolism in PA14 biofilm compared to transcriptome-free MiMICS.

From experimental data, plotted is the average percentage of total cells expressing a gene. From transcriptome-guided MiMICS and transcriptome-free MiMICS simulation data, plotted are the average percentage of total agents with active flux through the reaction encoded by a gene. Error bars represent one standard deviation from n = 7 experimental replicates and n = 50 simulation replicates. Experimental data provided by Dar and co-workers.

https://doi.org/10.1371/journal.pcbi.1012031.g004

As a result of the ABM rules enforced in transcriptome-guided MiMICS (Table 1), a heterogenous agent distribution simulating the denitrification and oxidative stress metabolic model state emerged within the spatially variable NO biofilm microenvironment (Fig 3). Consequently, reflective of the intracellular metabolic model states (Fig 2B), transcriptome-guided MiMICS predicted a relatively low number of cells (10% of cells) utilizing narG- and napA-encoded reactions compared to the cell population utilizing norB- and nosZ-encoded reactions (15% of cells) (Fig 4). This MiMICS prediction agreed with the experimental results reporting that relatively few cells expressed narG and napA, 1.5% and 4.3% of cells, respectively, compared to cell populations with high norB and nosZ expression, 12.1% and 11.8% of cells, respectively (Fig 4).

Interestingly, transcriptome-guided MiMICS did not accurately predict larger cell populations expressing nirS (20.4% of cells in experiment, 10% in simulation) relative to cell populations expressing narG, norB, and nosZ (specific percentages reported above) (Fig 4). While NO has been observed to also upregulate nirS expression³², hypothesized to bind NO required for NO transport to NorB for degradation [38], the current PA14 GENRE does not include this NirS NO-binding reaction. Thus, the NO-induced oxidative stress metabolic model did not predict active nirS-encoded reaction flux (Fig 2B), likely accounting for MiMIC’s prediction discrepancy of a larger nirS-expressing cell populations (Fig 4).

Only the transcriptome-guided MiMICS simulation predicted a NO biofilm microenvironment (Figs 3 and S5), which was hypothesized to induce oxidative stress in nearby cells. Accordingly, only transcriptome-guided MiMICS simulations predicted a small population of cells (0% in transcriptome-free simulation, 7.2% in transcriptome-guided simulation, and 12.5% in experiment) with flux through the katA-encoded oxidative stress reaction (Fig 4).

In terms of carbon metabolism, both transcriptome-free and transcriptome-guided MiMICS correctly predicted a majority of cells with flux through a succinate dehydrogenase reaction used in the tricarboxylic acid (TCA) cycle, encoded by sucC, and a glucose transport reaction encoded by gltB (Fig 4). The transcriptome-guided MiMICS prediction of a large population of cells expressing sucC (90% of cells) and gltB (84% of cells) is a result of the flux through sucC- and gltB-encoded reactions that was predicted to be highest in the aerobic metabolic model state (Fig 2A and 2B), which most agents represented because a majority of agents were likely exposed to high oxygen. Transcriptome-guided MiMICS accurately reduced cell populations (46.8% in transcriptome-free simulation, 15.0% in transcriptome-guided simulation, and 4.4% in experiment) with active flux through an acetate kinase fermentation reaction, encoded by ackA [22] (Fig 4). Related to energy metabolism, both MiMICS transcriptome scenarios predicted 100% of cells with active flux through the ATP synthase reaction, encoded by atpA, and CTP synthase, encoded by pyrG, compared to experimental populations expressing atpA (82.4% of cells) and pyrG (58.8% of cells) (Fig 4).

Potential inaccuracies in either the -omics data-integration method or the MiMICS model components rendered the transcriptome-guided MiMICS incapable of predicting the relatively low cell proportions utilizing gshB and aceA, which encode glutathione synthetase and isocitrate lyase, respectively, or cysA, which encodes a sulfate transport protein (Fig 4). In addition, neither MiMICS transcriptome version captured heterogeneous expression of genes related to virulence factor synthesis, which was expected as virulence factor synthesis likely does not contribute to the biomass synthesis GENRE objective function used during transcriptome data integration [37]. Overall, compared to the standard practice of using a transcriptome-free multi-scale framework, transcriptome-guided MiMICS improved prediction accuracy of the heterogeneous cell distributions with upregulated denitrification and oxidative stress metabolism.

To investigate the impact of the ability for cellular agents to dynamically adopt different metabolic model states, a MiMICS scenario was run with cells that remained fixed in a randomly initialized metabolic state, called ‘Fixed metabolic state MiMICS’ (S9 Fig). Fixed metabolic state MiMICS predicted homogeneous metabolic state niches localized to the initial location of the respective metabolic state and limited spatial mixing of different metabolic states. In contrast, the dynamic switching ability in MiMICS promoted microscale niches of cells with a heterogeneous distribution of metabolic states localized to the biofilm center. This result confirms that dynamic metabolic state adaptation to the metabolic microenvironment is essential to recapitulate the confined microscale niches of heterogenous metabolic states observed in the experimental biofilm.

MiMICS captured spatial relationships of intracellular metabolism in PA14 biofilm

To quantitatively evaluate the validity of transcriptome-guided MiMICS predictions in space, a bulk neighborhood spatial correlation analysis between genes was performed, similar to the analysis of experimental data by Dar and co-workers [5] (Fig 5B). Genes were compared which had both positive experimental expression values and simulated active flux values of reactions encoded by the respective gene. Shown in Fig 5B, MiMICS accurately predicted the spatial correlation for 33 gene pairs, and incorrectly predicted the spatial correlation for 22 gene pairs. Specifically, MiMICS accurately predicted that denitrification genes (i.e. narG, nirS, norB, nosZ) were positively correlated with one another, and positively correlated with the oxidative stress gene katA (Fig 5A and 5B). Conversely, MiMICS predicted napA, which encodes nitrate reductase, was positively correlated with the remaining denitrification genes narG, nirS, norB, and nosZ. However, this correlation was not observed in the experiment. Upon closer inspection, MiMICS predicted patches in the biofilm where neighboring cells utilized napA-, nirS-, norB-, and nosZ-encoded reactions (S6 Fig). In contrast, cells expressing napA in the experiment were more dispersed compared to more spatially confined biofilm regions of cells expressing nirS, norB, and nosZ (S6 Fig). MiMICS prediction discrepancies in the napA spatial correlations suggest alternate mechanisms exclusively modulate napA expression, such as extracellular phenazine secretion [39]. Indeed, phenazine synthesis genes phzE1 and phzM were expressed in the ten-hour PA14 biofilm (Fig 4). In addition, MiMICS accurately predicted denitrification and oxidative stress genes (i.e. narG, nirS, norB, nosZ, katA) were anticorrelated with the carbon metabolism genes gltB and sucC (Fig 5A and 5B). One aspect MiMICS did not accurately predict was the spatial correlation of atpA, encoding ATP synthase, with all other genes (Fig 5B), which motivates improvements in future transcriptome-guided MiMICS simulations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. MiMICS captured spatial relationships of intracellular metabolism in PA14 biofilm.

(A) Representative xy projections of PA14 biofilms from experiments and transcriptome-guided MiMICS simulations. Cells plotted are located near the z = 0 μm surface. Colored cells have high expression of the gene listed (experiment) or high reaction flux encoded by the gene listed (simulation). Circled areas highlight regions of interest where there is an anticorrelation between TCA cycle metabolism with denitrification and oxidative stress metabolism. Scale bar represents 20 μm. (B) Neighborhood gene spatial correlation comparison between experiment and transcriptome-guided MiMICS simulation. Spatial correlation between gene pairs was assessed by a Pearson correlation, where +1 and -1 value correspond to a strong positive and strong negative spatial correlation, respectively. The experimental and simulation Pearson correlation values are plotted in the upper right and lower left triangle of each square, respectively. Simulation Pearson correlation values were determined from 50 simulation replicates. Experimental data was reconstructed from Dar and co-workers.

https://doi.org/10.1371/journal.pcbi.1012031.g005

MiMICS revealed functional dynamics of denitrification and oxidative stress metabolism in PA14 biofilm

This spatial transcriptomics experiment required biofilm fixation, preventing in-situ dynamic monitoring of gene expression within the biofilm sample. Thus, this spatial transcriptomic method was unable to characterize the temporal events that underline the PA14 biofilm gene expression distributions and spatial patterns. In contrast, transcriptome-guided MiMICS dynamic simulations were used to monitor and quantify spatiotemporal shifts of PA14 biofilm metabolism connected to the metabolite microenvironment.

Agents were randomly initialized at t = 0hrs, nucleating the simulated PA14 biofilm which was grown over a period of ten hours (Fig 6). After approximately nine hours of simulated PA14 biofilm growth, agents exposed to the emergent microaerobic biofilm environment (i.e. local oxygen below the [O₂]_t value) switched from the aerobic to denitrification metabolic model state (Fig 6A and 6B). In accordance with the aerobic and denitrification metabolic model states (Fig 2A), cell populations with active flux through narG-, nirS-, norB-, and nosZ-encoded denitrification reactions increased over time (Fig 6A and 6B). The stochastically-chosen agents that executed the norB-limiting denitrification metabolic model (Table 1) secreted NO, generating the NO biofilm microenvironment (Fig 6A). The maximum level of NO, ~4 μM, below toxic levels [23], was localized in the biofilm center (Fig 6A). At ten hours, agents that sensed local extracellular NO above the [NO]_t value switched to the oxidative stress metabolic model state, resulting in increased cell populations with flux in the oxidative stress reaction encoded by katA (Fig 6).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. MiMICS revealed the microscale dynamics of denitrification and oxidative stress intracellular metabolism, and extracellular nitric oxide (NO) in the PA14 biofilm.

(A) Representative xy projections at a z = 0 μm slice of a simulated PA14 biofilm colony developing over time. Agents are colored according to their assigned metabolic model state. Agents are also colored according to active flux through reactions encoded by denitrification genes (narG, nirS, norB, and nosZ) and encoded by the oxidative stress gene katA. Plotted are the corresponding xy projections of oxygen and NO concentrations at the z = 0 μm slice. Scale bar represents 10 μm. (B) Quantified dynamics of the average number of agent’s assigned a metabolic model state, and the average number of agents with active reactions flux encoded by denitrification and oxidative stress genes. Dynamics of the average extracellular oxygen and NO in the biofilm are reported. Error bars represent one standard deviation from 50 simulation replicates.

https://doi.org/10.1371/journal.pcbi.1012031.g006

Due to the ABM rules that were enforced (Table 1), a heterogeneous distribution of cells using the denitrification and oxidative stress metabolic models emerged within microaerobic and spatially variable NO biofilm microenvironments. Only the denitrification +/-NO metabolic model states predicted active flux through narG- and nirS-encoded intracellular reactions, which are upstream of NO synthesis and required for NO secretion (Figs 2B and 6A). Conversely, denitrification and oxidative stress metabolic model states predicted active flux through norB- and nosZ-encoded intracellular reactions, which are downstream of NO synthesis and used for NO degradation (Figs 2B and 6A). As a result of these intracellular differences among metabolic model states, fewer cells in the biofilm were predicted to use narG- and nirS-encoded reactions compared to cells using norB- and nosZ-reactions (Fig 6B).

Relating these MiMICS dynamic predictions to possible biological functions, previous studies similarly observed cells downregulated narG expression and upregulated expression of nirS and norB, the latter encoding a NO degradation reductase, in presence of NO [40,41]. This NO-induced regulation of different genes in the denitrification pathway was hypothesized to arrest NO synthesis, promote NO degradation, and maintain extracellular NO within non-toxic concentrations [40,41]. Functioning similarly, MiMICS predicted NO signaled agents to reduce flux through narG and nirS-encoded NO synthesis reactions, which decreased NO synthesis and secretion (Fig 6). Due the decreased NO secretion, and extracellular NO gradients dissipated at 9.9 hours in the biofilm colony shown in Fig 6A. In contrast, possibly to promote degradation of extracellular NO, agent populations with active flux through norB- and nosZ-encoded NO degradation reactions continued to increase (Fig 6).

Next, possibly due to the reduced NO concentrations and expanding microaerobic biofilm regions, agents switched to a denitrification +/-NO secretion metabolic model state, again promoting a NO-rich biofilm microenvironment at 10 hours in the biofilm colony in Fig 6A. Thus, MiMICS predicted an oscillating, non-toxic NO biofilm microenvironment (Fig 6). This prediction is supported by previous studies that have measured temporal oscillations in extracellular NO in low oxygen P. aeruginosa cultures performing denitrification [42,43]. Furthermore, temporal oscillations in expression of mRNA encoding denitrification reductases has been observed in E. coli [44]. In MiMICS, the predicted oscillations in extracellular NO were a result of cellular agents differentially regulating intracellular reactions in the denitrification pathway in response to extracellular oxygen and NO, which may function to maintain a non-toxic NO biofilm microenvironment. Future simulations can systematically perturb genes and parameters in MiMICS to test their effect on predicted outcomes, such as NO metabolism in the biofilm. For example, in silico gene knock-outs, such as ΔnirS that increases NO secretion, or varied extracellular nutrient concentrations may be screened for increased NO-induced P. aeruginosa cell death as a potential therapeutic strategy. Simulation readouts, such as biofilm structures, gene expression, or number of live/dead cells, can be quantified and compared to the experiment.

MiMICS computational performance

Performing GENRE simulations for individual cellular agents is advantageous to predict emergent metabolic heterogeneity within a multicellular system but can be computationally expensive. Execution of one MiMICS simulation time step on one central processing unit (CPU) resulted in a runtime on the order of 10 minutes for 10,000 cellular agents. To improve this MiMICS runtime, GENRE simulations for each individual agent were divided across multiple CPUs using parallel computing, implemented with Python Multiprocessing. This parallel computing strategy reduced MiMICS runtime by an order of magnitude for 10,000 cellular agents (S7 Fig). This computational performance is similar to simulation runtimes reported by the BacArena framework [17] for 1000 agents, where MiMICS is the first framework to report efficient runtimes for up to 10,000 single-cell agents. The number of metabolic model states input into MiMICS did not considerably impact MiMICS runtime (S7 Fig).

Genome-scale network reconstruction (GENRE)

Pseudomonas aeruginosa PA14 metabolism was simulated with a published P. aeruginosa PA14 genome-scale metabolic network reconstruction (GENRE) (iPau21) [12,37]. All GENRE simulations were performed in Synthetic Cystic Fibrosis Sputum Medium (SCFM). A few reactions were updated in the PA14 GENRE to represent the denitrification, fermentation, nitric oxide secretion, and oxidative stress metabolic processes observed in experiments. PA14 GENRE simulations in anaerobic + nitrate were initially infeasible (S8 Fig), which did not reflect expected denitrification biomass synthesis processes [22]. In anaerobic SCFM conditions, ubiquinol-9 (UQ₉), a preferred respiration cofactor required in the PA14 GENRE biomass objective function [37], synthesis was found to be blocked which prevented biomass growth (S8 Fig). Three oxygen-independent hydroxylation reactions of UQ₉ precursor metabolites were added with H₂O as the reactant in the place of oxygen, a hypothesized oxygen-independent reaction [46,47,48]. Oxygen-dependent UQ₉ synthesis reactions were set to be irreversible to prevent the generation of intracellular oxygen. To promote feasible L-arginine fermentation, reactions encoded by arcB and argH were set as reversible and irreversible, respectively, in agreement with the expected arginine deiminase pathway [49]. These updates to the PA14 GENRE promoted feasible biomass growth rates and flux through the respective reactions in anaerobic conditions (S8 Fig). Exchange and transport reactions were added for the denitrification intermediates nitric oxide (NO) and nitrous oxide (N₂O). NO toxicity was represented as a sink reaction for NO and iron [23] and an oxidative stress KatA reaction was represented as a sink reaction for NO [23]. Of note, when cells were exposed to an extracellular NO microenvironment in MiMICS, NO uptake was enforced by fixing the GENRE upper bound for NO exchange at a negative flux value, which was calculated from the cell’s local extracellular NO concentration (Eq 1).

In MiMICS, each agent’s PA14 metabolic model was optimized for a biomass synthesis objective function using flux-balance analysis (FBA) [9]. Each agent’s metabolic model exchange reaction bounds were constrained by the agent’s local oxygen, nitrate, glucose, and nitric oxide concentrations. The agent’s local patch metabolite concentrations and metabolic model lower bound of exchange reaction flux were converted with the following equation: (1) Where f_M is the metabolite flux that set the metabolic model exchange flux (mmol/(gdWt *hr), c_M is the agent’s local patch metabolite concentration (mM), v_patch is the patch volume unoccupied by the bacteria (L), dt_ex is the metabolite uptake time step (estimated as 0.05s in Table 2), and b is the agent biomass (g). v_patch was assumed as 1e^-16 L.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. MiMICS parameters and estimated time scales.

https://doi.org/10.1371/journal.pcbi.1012031.t002

Construction of transcriptome-guided PA14 metabolic models

The PA14 GENRE was integrated with a published spatial transcriptomics PA14 biofilm dataset [5] using the established RIPTiDe algorithm [11]. The dataset consisted of the mRNA expression of ~292,000 PA14 biofilm cells grown in SCFM, and fixed at either 10hr (n = 7) or 35hr (n = 3) growth time points [5]. For data integration into the PA14 GENRE, gene expression values were only used for the 47 of the 105 genes measured in the experiment that also encoded reactions in the PA14 GENRE (S1 Table). For data integration into the PA14 GENRE, lower bounds for exchange reactions were constrained according to SCFM concentrations (S2 Table), converted using Eq 1. In addition, upper bounds for exchange reactions were set to +1000 to allow for unconstrained production of exchange metabolites.

A UMAP Leiden clustering analysis (scanpy v.1.7.0), data preparation and parameters described elsewhere [5], identified distinct gene expression cellular states within the spatial transcriptomic PA14 biofilm dataset (S1 Fig). Metabolic states for each cluster were assigned based on the highest ranked genes within the cluster. Median gene expression was extracted for the top 9 clusters, which captured 91% of biofilm cells, that were integrated into the PA14 GENRE using RIPTiDe [11] to generate 9 UMAP-identified PA14 metabolic model states. A fractional growth rate of 0.7 was used in RIPTiDe. Reactions not shared by all UMAP-identified metabolic model states, termed non-consensus reactions, and predicted exchange metabolite fluxes were compared and grouped based on similarity to identify unique 10-hour and 35-five PA14 metabolic state models (S1 Fig). Due to the wide diversity of metabolic states between the 10hr and 35hr time points, the initial scope of this study was on metabolic states emerging within PA14 biofilms grown for 10 hours. This UMAP-informed method generated the 10hr biofilm aerobic PA14 metabolic model state that was input into MiMICS.

A manual selection method was used to extract gene expression values of denitrification and oxidative stress metabolic subpopulations that were not identified in the UMAP Leiden clustering analysis. The ~22,000 denitrification cells classified in the clustering analysis were categorized into having high expression of narG and/or napA (S2 Fig). The narG- and napA-expressing denitrification-classified cells were further classified into possessing a non-limiting or limiting denitrification pathway (S2 Fig). Cells with a non-limiting denitrification pathway did not possess a rate limiting denitrification transcriptional step (i.e. expression of narG < = nirS < = norB < = nosZ). Cells categorized with limited denitrification metabolism had expression of a gene encoding a denitrification reductase that was larger than the expression value of the gene encoding the subsequent denitrification reductase gene. From cells classified with non-limiting and limiting denitrification metabolism, the median gene expression was extracted and integrated into the PA14 GENRE using RIPTiDe. The fraction of the optimal biomass was set as 0.55 used in RIPTiDe for these denitrification subpopulations [30]. The method generated four PA14 denitrification metabolic model states: no limiting denitrification step, nirS-limiting, norB-limiting, and nosZ-limiting. Predicted exchange fluxes of each of these denitrification limiting metabolic model states are shown in S2 Fig. This method generated the denitrification +/- NO secretion metabolic model states that were input into MiMICS.

An oxidative stress PA14 metabolic model state was generated by integrating the PA14 GENRE with median gene expression values extracted from the ~22,000 denitrification-classified cells with high katA expression. To simulate the exposure of an NO-induced oxidative stress environment, during transcriptome integration into the PA14 GENRE, the upper bound for the NO exchange reaction was set to a negative NO flux value (converted from 20 μM NO) to enforce NO uptake. The fractional optimal biomass was 0.7 used in RIPTiDe. This method generated the oxidative stress PA14 metabolic model that was input into MiMICS.

Metabolite reaction-diffusion model

Metabolite concentrations were simulated with the built-in HAL partial differential equation (PDE) alternating direction implicit (ADI) method solver. In HAL, the metabolite grids were defined as 3-dimensional 230 μm x 230 μm x 40 μm grids divided into 2 μm x 2 μm x 2 μm patches (Table 2). Metabolite grids were constructed for extracellular oxygen, nitrate, nitric oxide, and glucose.

Metabolite concentrations were initialized in each patch based on a Synthetic Cystic Fibrosis Sputum Medium (SCFM) recipe (S2 Table ) [29]. The initial condition for the metabolite PDEs was assumed to be spatially uniform: where c_initial was the initial metabolite concentration (mM).

Before metabolite diffusion simulations, a metabolite reaction term was discretely applied at each agent’s respective grid location, representing the local metabolite consumption or production by that agent. The reaction term consisted of first converting the agent’s metabolic flux, predicted by the agent’s metabolic model, to a concentration (Eq 1), which was then passed to update the metabolite concentration at the agent’s grid location (Eq 2). (2) where c_f is the updated metabolite concentration (mM) at the agent’s location, c₀ is the initial metabolite concentration (mM) at the agent’s location, and Δc is the agent’s metabolic model predicted metabolic flux converted to a concentration (mM). With a consideration towards the thousands of agents, each with multiple predicted metabolic reaction rates, this discretely applied metabolite reaction rate approach was used to improve computational efficiency.

Next, diffusion was simulated using HAL’s ADI diffusion equation solver method. The governing equation for metabolite diffusion was as follows: (3) where c is the metabolite concentration (mM), t is time (s), x is the x-dimension, y is y-dimension, z is the z-dimension, and D is the metabolite diffusion coefficient (cm²/s).

No flux boundary conditions were applied in the z-dimension: where z = 0 was at the glass slide and z = L was at the air-aqueous interface.

Periodic boundary conditions in the x- and y-dimensions were applied. The biofilm region was defined as R_biofilm. Regions outside of R_biofilm assumed constant metabolite concentrations. This assumption was attributed to higher metabolite diffusion coefficients occurring in the aqueous phase compared to metabolite diffusion coefficients the biofilm. In addition, the aqueous nutrient-rich media was refreshed every four hours in the experiment [5], so constant nutrient concentrations in the aqueous phase were assumed. To apply this assumption, for every metabolite diffusion step, the metabolite concentrations were set to the initial concentration in regions outside of the biofilm:

Aqueous diffusion coefficients were multiplied by 0.6 for light gases and 0.25 for organic solutes to calculate the diffusion coefficient in a biofilm (Table 3) [53]. Biofilm diffusion coefficients were scaled according the recommended HAL equation for PDE stable solutions [24]: (4) where D is the metabolite diffusion coefficient (cm²/s), dt_D is the diffusion time step (seconds), and dl is the patch length (cm) (2 μm in Table 2). Gaseous and carbon metabolite diffusion coefficients were calculated according to the respective metabolite diffusion time scales, dt_{D_gas} and dt_{D_carbon}, in Table 2.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Metabolite diffusion coefficient values used to solve metabolite PDEs.

https://doi.org/10.1371/journal.pcbi.1012031.t003

Agent-based model

The agent-based model (ABM) was constructed in the established Java-based platform Hybrid Automata Library (HAL) [24]. The ABM world was defined as a 3-dimensional 230 μm x 230 μm x 40 μm grid divided into 2 μm x 2 μm x 2 μm patches. Periodic boundary conditions were set in the x and y-directions. An individual agent represented an 8 μm³ single-cell P. aeruginosa bacterium. The agent class was built using HAL’s off-lattice SphericalAgent3D agent class. Each agent was assigned a unique index number. Bacteria agents were randomly initialized at z = 0 μm with a random biomass in the range of 1e^-12 to 2e^-12 g and a random directional angle. Biomass growth was calculated with an exponential growth rate law [5]: (5)

Where b is the updated agent biomass in grams, b₀ is the initial agent biomass in grams, μ is the growth rate informed from the agent’s optimized metabolic model (hr^-1) and dt_growth is the growth time step (hr) (Table 2). If an agent biomass grew above the maximum biomass threshold (Table 4), cell division occurred, and a new daughter agent was placed in a neighboring patch of the mother agent. The mother agent biomass was randomly divided between the mother and daughter agent. The daughter agent’s directional angle was set to within 10 degrees of the mother agent’s angle. The daughter agent’s metabolic state was set as the mother agent’s metabolic state. An agent became inactive, representing cell death, in the agent-based model simulation if the predicted growth rate was 0.0 hr^-1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Parameter values used for bacteria agents.

https://doi.org/10.1371/journal.pcbi.1012031.t004

Simple ABM rules for P. aeruginosa surface motility were implemented based off previous observations of physical cell-cell interactions and collective swarming movement of piliated P. aeruginosa [36,55]. Pili surface motility speed was previously reported as 10 nm/s [51]. Thus, for a 5 minute simulation time step, one piliated bacteria agent was estimated to move 3 μm or 1.5 patches. For each five-minute simulation time step, the surface motility ABM function was performed once for each piliated bacteria agent. To determine the agent’s direction of movement, an agent on the surface (z = 0 μm) with at least one unoccupied and occupied neighbor patch set its orientation angle equivalent to a random neighbor agent. Agents moved to the new location if there were fewer than 20 cells in the surrounding region of the new location. For motility simulations, all bacteria agents in the ten-hour biofilm growth period were assumed to express pilA (S1 Fig) and perform the surface motility ABM function. For non-motile simulations, the motility ABM function was removed for all agents. In addition, to prevent agent overlap, the built-in HAL cell mechanics algorithm was executed for each agent using force scalar and friction coefficient parameters listed in Table 3.

Each agent was assigned their own metabolic model to simulate metabolism. Four PA14 metabolic model states were available for agents to choose from, specifically, the aerobic, denitrification, denitrification +NO secretion, and oxidative stress metabolic model states. Agents decided their metabolic model state by comparing their local oxygen and NO concentrations to the respective metabolite threshold values, [O₂]_t and [NO]_t (Tables 1 and 4). For agents exposed to low oxygen and low NO environments, a stochastic parameter, R_n, generated a number between 0 and 1. These cells decided between a denitrification with or without NO secretion metabolic model state by comparing their R_n value to a threshold value, 0.06, which was informed from the experimental percentage of denitrification cells with a norB-limiting expression profile. After agents checked these conditions, the agents were given a metabolic state assignment in the ABM, represented by an integer value. Used to assign the correct metabolic model state to each agent in Python, this metabolic model state integer value corresponded to the order in a list of metabolic model states input into MiMICS.

A biologically-relevant 10 minute time delay, dt_mRNA, for agent’s to switch to a new metabolic model state upon sensing a new environmental cue was implemented, corresponding to the 10 minute delay for synthesis and detection of mRNA transcription upon a cell sensing a new environmental metabolite cue [52] (Table 2). Thus, when an agent encountered a new metabolite cue that induced a metabolic model state switch, the time the bacteria agent occupied that new metabolite environment was recorded. Once this recorded time exceeded 10 minutes [52], the bacteria agent was assigned the new metabolic model state. For dividing agents, this recorded time was equally split between the mother and daughter agents, representing observations of detected mRNA amounts in dividing cells [52].

Parameterization

The oxygen threshold parameter, [O₂]_t, was varied from 0.19–0.23 mM. The mean absolute error between experiment and MiMICS outputs was calculated from 21 replicate simulations performed for each [O₂]_t parameter value. The total MiMICS model error was calculated from the summation of the error from MiMICS predictions between two experimental outputs: the percentage of denitrification classified cells, and the percentage of oxidative stress classified cells. Denitrification cells were classified to have expression of one or more the following denitrification genes: nirS, norB, nosZ, and no expression of the oxidative stress katA gene. Oxidative stress cells were classified to have expression of the oxidative stress katA gene. The [O₂]_t parameter value with the smallest error relative to the experiment was used for subsequent simulations.

MiMICS specifications and code availability

In MiMICS, a Py4J JavaGateway Server was implemented to interface the Java-based HAL platform and Python-based metabolic model. Through this Py4J JavaGateway, functions are called from Python to initialize and run the ABM. Additionally, Python functions called through the JavaGateway retrieve and pass information between Python and Java and save spatial agent and metabolite information at each simulation time point. MiMICS was initially developed on a Mac computer (macOS Catalina v.10.15.7, 1.1 GHz Dual-Core Intel Core i3) using IntelliJ IDEA and 1–2 central processing units (CPUs). All MiMICS simulations were executed on Rivanna High-Performance Computing (HPC) System at the University of Virginia. MiMICS Java files were compiled as a JAR file to run on the Rivanna HPC System. To reduce simulation runtime, MiMICS split metabolic model calculations for each agent across 35 CPUs using Python multiprocessing v.0.70.14. One MiMICS simulation of the ten-hour biofilm growth period on Rivanna HPC System across 35 CPUs took 40 minutes and required 40 GB of memory. When appropriate, multiple MiMICS simulations (e.g. for simulation replicates, parameter sweeps) were executed across multiple computing nodes to prevent multiple JavaGateway servers running on a single computing node. MiMICS simulations used Python v.3.9, cobra v.0.29.0, py4j v.0.10.9.7, Java v.1.8.0, HAL v.1.1.0. MiMICS source code and detailed user guide can be found at: https://github.com/tracykuper/mimics.

Statistical analysis

Statistical analysis was performed with the scipy stats python package (v.1.10.0). Variation among biofilm structure metrics among experimental and simulation conditions was performed with analysis of variance. Pairwise comparisons of individual group means were performed using a Tukey post hoc analysis. Values of p < 0.05 were considered statistically significant.

Spatial neighborhood analysis

Bulk neighborhood analysis of simulation data was performed similarly to the experimental neighborhood analysis [5]. For each gene of interest, agents with active flux encoded by the reaction in the 99^th percentile were selected as “center” agents. Using the 3D centroid agent coordinates, distances between center agents and nearby agents within 3 μm were calculated. The closest five agents to the center agent were selected as neighbor agents. Neighbor agents were grouped together and the average reaction flux of the neighbors was calculated. This average neighbor reaction flux was divided by the corresponding population average reaction flux, omitting the center agents in the population. For each pair of genes, a Pearson correlation was calculated to determine spatially correlated genes among 50 simulation replicates.