Bacteria Strain and Growth Condition

Bacterial Strain: We utilized the dual labeled Pseudomonas aeruginosa, with a QS reporter. The strain, PAO1(wt)::lasB-gfp(ASV)::Plac-lasR::rfp, was generously provided by Michael Givskok23. Furanone C-30 was acquired from Sigma-Aldrich Co., and used at a final concentration of 100 μM.

P. aeruginosa biofilms: Biofilms were set up on MatTek plates, at 37 °C using Columbia browth. Furanone was added at time of 24 hour after biofilm development. Biofilm images were collected on a Zeis510 Meta Conforcor3 Laser Scanning Microscope at 25, 27, and 28.5 hours post-seeding.

Bacterial Image Processing

We detect the intensity of each bacterial image acquired from the confocal laser. Next, we compute the number of bacteria that express rfp and gfp, respectively. We then compute the relative concentration of the QS signal by using the following formula:

where N rfp and N gfp represent the number of bacteria that express rfp and gfp, respectively.

Definition of the Bacterial Genotypes and Regulons

Bacterial genotypes: Specifically, we make use of three naturally occurring genotypes in regards to QS genes 39,40 (i) functional QS systems sensitivity to QSI ( QS + ), (ii) deletion of the QS machinery ( QS − ), (iii) or modified QS systems that are functional but resistant to QSIs (QSI resistant). QS − cells that have access to the benefits of QS, but do not incur the metabolic costs of production, should be at an adaptive advantage relative to their QS + counterparts and can be characterized as “cheaters” within the community. We assume that the strains do not differ in other regions, and consequently there are no fitness advantage incurred by variability in other genomic regions.

Bacterial regulons: Signaling via QS can lead to the production of multiple types of goods (Fig. 3). We separate these into four categories: (i) non-beneficial goods that provide no fitness advantage to producers or neighboring cells; (ii) private goods that provide a fitness advantage to producers but not neighboring, (iii) quasi-public goods that provide a fitness advantage to producers and neighboring cells but remain exclusively accessible to their producers, and (iv) public goods that provide a fitness advantage to producers as well as neighboring cells. The competitive index between QS+ and QS− depends the production costs and benefits of QS, and the availability of QS products to QS− cells.

Bacteria Population Model

The intracellular biochemical pathways giving rise to quorum sensing have been well studied from both molecular and systems biological perspectives41,42,43,44. Yet, the endogenous mechanisms of such signaling systems are inextricably bound to exogenous processes, cues, and constraints suggesting that, in isolation, intracellular models are unrealistic45,46. Consequently, recent efforts have been put forth to contextualize the intracellular signaling mechanisms within the scope of intercellular interactions47,48,49,50,51,52,53,54,55. Nevertheless, this line of work has, hitherto, been limited to interactions among small groups of cells and has neither considered characterizations of molecular communication at a higher level of abstraction, nor differentiated cell types capable of behavioral variances including cooperation, exploitation, and pathogenesis.

To this end, we present a new computational model for bacterial community dynamics capable of capturing highly emergent behaviors found at the population level. We draw from the paradigm of agent-based modeling whereby each cell is outfitted with its own set of equations governing metabolic and communication processes56,57,58,59,60. Our model domain consists of a 3-dimensional micro-fluidic environment governed by the laws of diffusion and volume exclusion. The model incorporates molecular signaling, nutrient limited metabolism, exoproduct formation, and virulence processes that are crucial in the development and evolution of mixed biofilm communities (see Supplementary Table 2 for model parameters).

In our model, we assume the QS regulatory network of Pseudomonas aeruginosa has two feedback loops (see Supplementary Fig. 9). The LasR − AHL complex up-regulates the expression of both lasR and lasI genes, generating even more signal molecules (AHL) and receptors (LasR), thereby forming a positive feedback loop. To model the QS system, we use a set of ordinary differential equations proposed in refs 47 and 61, and then extend the ODE system to incorporate the effects of QSI. Particularly, we focus on three components of the QS system that are potential targets for QSI therapy: (1) LasI signal synthase (i.e., signal generation) (2) Extracellular AHL molecules (i.e., signal accumulation) (3) LasR antagonist (i.e., signal reception)32,33.

Furthermore, we simulate bacterial growth kinetics based on the classic Monod model62 by considering multiple nutrient sources. The production of all inter or intracellular proteins are all constrained by nutrition availability in the ambient microenvironment. Nutrients and drug molecules diffuse across the space following Fick’s low while using different diffusion coefficient for different layers (see Fig. 1). Physical collisions are handled using CUDA platform which is capable of calculating the overlapping area of thousands of particle pairs (i.e., cells and EPS) and shoving them simultaneously.

Quorum Sensing Model

The regulatory network of the LasR/R quorum sensing system has two feedback loops. Based on the ODE-models proposed in refs 47 and 61, we have the following equations for the luxIR QS system:

where [X] denotes the concentration of a particular molecular species X. In our formulation, A stands for AHL, R is LasR homologs, RA is the LasR − AHL complex and C is the dimerized complex. c A and c R account for the basal level transcription of A and R, respectively. The values and references of the parameters are adapted from a general LasIR system47.

To give some intuition, the first term c A of Eq.(2) describes the basal level transcription, the second term captures the positive feedback loop regulated by the dimerized complex C and the third and forth terms describe the AHL concentration changes caused by the binding and unbinding reactions of AHL and LasR receptor, respectively. The last term describes the degradation of AHL. Eqs (4) and (5) describe the binding reaction of AHL and LasR receptor, as well as the dimerization process of the binding product [RA], respectively.

Quorum Sensing Inhibition Model

Recently, a large number of chemical compounds have been screened and many of them show QS inhibition effects targeting various components of the QS system. In this paper, we focus on three components of the QS system that are potential targets for QSI therapy: (1) LasI Signal Synthase (i.e., signal generation) (2) Extracellular AHL molecules (i.e., signal generation) (3) LasR antagonist (i.e., signal reception)32,33.

(1) Inhibition of LasI signal synthase. AHL synthase catalyses the formation of an amide bond between the homoserine lactone ring from S-adenosylmethionine (SAM) and the acyl chain from the acyl acyl-carrier-protein (acyl-ACP). After catalysed, the product of AHL and by-products of holo-ACP and 59-methylthioadenosine (MTA) are generated and released63. Many chemical compounds have been shown to be effective reducing LasI activity14. Recently, a new compound trans-cinnamaldehyde, which has no chemical structure similarity to AHLs or AHL analogues, is shown to inhibit AHL synthases by molecular docking studies34. To have a general form in our model, we modify eq. (2) to add the lasI inhibition; also, a new equation for the drug degradation is also added:

where [QSILasI] stands for the concentration of QSI molecules that inhibits luxI expression, k 6 is its degradation rate.

(2) Extracellular AHL: AHL molecules can diffuse in and out of the bacterial cell freely. Therefore, once they appear in the extracellular environment, they are potential targets for destruction or inactivation. The AHL-degrading enzyme, AHL-lactonase (such as AiiA or AiiM), has been reported to deactivate the bacterial virulence through hydrolysis of the lactone ring of AHL64. Accordingly, eq. (6) needs to be modified to describe the AHL-degradation by AHL-degrading enzyme; similarly, we also need to add a new equation describing the drug degradation:

where [QSI d ] stands for the concentration of QSI which degrades AHL. In these equations, k cat represents the maximum rate achieved by the system, at maximum (saturating) substrate concentrations. The Michaelis constant K M is the substrate concentration at which the reaction rate is half of k cat , and k 7 is the drug degradation rate.

(3) AHL signal reception: The most promising mechanism for inhibiting LasR activation is achieved through the use of AHL analogues that act as antagonists for the native AHL (i.e., 3O-C12-HSL for LasIR system). These molecules are likely to be similar in structure to the natural AHL and compete for LasR-receptors binding (some can bind covalently). Accordingly, eq. (3) needs to be modified to:

where [QSI anlg ] stands for QSI which inhibits the LasR activation, [RQSI anlg ] is the binding product of LasR and [QSI anlg ]. Also, we need to add two more equations to describe the dynamics of [QSI anlg ]:

where k 10 is the degradation rate of QSI anlg . We assume that the AHL analogues have the twice the affinity as the native AHL to the LasR receptor, therefore, the binding reaction rates satisfy k 8 = 2k 1 .

Cell Growth

Monod introduced the concept of single nutrient controlled kinetics to describe microbial growth62. The kinetic relates the specific growth rate (μ X ) of a bacterium cell mass (X) to the substrate concentration (S). The kinetic parameters, maximum specific growth rate (k X ) and substrate affinity (K S ), are assumed to be constant and dependent on strain, medium, and growth conditions (e.g. temperature, pH). Also, we consider a second nutrient source, and add a new term Q. On the other hand, when cells are metabolically active, but not growing or dividing, they may still take up substrate. To address this, a maintenance rate (m) is generally used to describe the reduction, and Monod’s model can be improved as follows:

EPS Production

We model the production of EPS as a function of the intracellular dimerized LasR − AHL complex C and incurring a cost on the carbon substrate (S) in Eq. (20):

where k EPS is the maximum EPS production rate, and represents nutrition limitation. Using this model, we expect a maximum rate of EPS production once cells reach a quorum through communication (i.e. the level of AHL in the extracellular environment passes a given threshold).

Digestive Enzyme (DGE) Production

We present a simplified and generic model of a diffusible public good motivated the production and utilization of the iron-chelating siderophore, pyoverdine. The production of DGE is also regulated by the intracellular dimerized LasR − AHL complex (C) a maximum production rate (k DGE ),

and the amount of complex Q (e.g., chelated iron) available in the extracellular environment is

where β is the chelation rate, DGE in and DGE ex are the intracellular and extracellular DGE concentration, respectively.

Metabolic Burden and Nutrition Consumption

Also, we need to take into account the cost of generating QS-related molecules, EPS, and digestive enzymes. Specially, to balance the limiting nutrition, we introduce utilization coefficients U X , U DGE , U EPS , U A , U R to model the metabolic burden and consumption of substrate nutrition to produce general cell mass, digestive enzymes, EPS, AHL molecules, and QS receptors, respectively:

Therefore, higher cell densities can lead to a decreased growth rate as well as production rate of various cellular products in a nutrition-limited environment.

Mass Volume, Radius, and Physical Interactions

Similar to60, in our simulation platform, a bacterial cell structure is compartmentalized with inner “biomass” core consisting of all intracellular material, and an outer layer consists of the capsular EPS. For each compartment, a density ρ i , the mass m j , and volume V j is updated at each simulation step according to the following equation:

Then, the radii of the inner biomass and entire cell are calculated using the volumes V j and V jtotal respectively. When the volume of a cell grow to be twice of the regular size (division threshold), it divides into two daughter cells. Similarly, when the capsular mass is above a “separate threshold”, the capsular mass fall off and form a new particle, such as the EPS.

Besides chemotaxis, cell growth, division, and shrinking are all the sources of movements. We use a physical collision kernel powered by CUDA technology to resolve any pair-wise agent overlap.

Cell-Cell Interaction Network and QS Signaling Molecules

Our network model considers the intracellular molecular information of the QS system, as well as the extracellular physical diffusion limit.

To replicate cell-cell interactions, in our model, a directed link from bacterium A to bacterium B is established under two conditions (see Fig. 8):

Figure 8 Bacteria QS-based collaboration network. A directed link is formed from cell A to cell B if cell B lies within the influence range, D if , of cell A and both cells have activated QS systems. The influence range, D if , of a cell depends largely on the intracellular concentration of the QS activity regulator complex, C. This is demonstrated by the bacterial color scale which represents intracellular [C]. The arrow from A to B indicates that A can influence B; the blue circles represent the radius of QS influence. Full size image

(I) Bacterium B must be within a diffusion-limited signal influence range D if of bacterium A.

(II) Bacteria A and B must be QS active at the same time.

The first condition accounts for the spatially constrained nature of molecular diffusion. The signaling molecules produced and secreted by bacteria have a strong impact only within a range, D if , depending on production rate, diffusivity, and decay rate. The second condition ensures that the QS system of both bacteria is active; this is realized by comparing against a concentration threshold of the intracellular QS regulator. To account for different levels of signal production of spreader bacteria, D if is defined as:

where α is a scaling factor for the overall influential distance, β is the saturation factor for the intracellular QS activity regulator complex C, and is the ratio of the effective diffusion (D) coefficient in the complex extracellular environment normalized to the free diffusion case (D 0 ). It is important to emphasize that the influence range, D if , differs across cells as a function of intracellular QS activity regulator complex, C. Fig. 8 illustrates an instance of a collaboration network formed across a small number of cells.

Physical Metrics to Quantity Biofilm Structure

Biofilm average thickness: The average biofilm thickness is measured based on the distance between the substrate and the outermost cell over the space.

Biofilm roughness: Biofilm surface roughness is defined by the following equation:

where L f is the average thickness, L fi is the ith value of the thickness matrix in the ith grid.

Metrics to Quantify Interaction Dynamics and Parameters

Once the cell-cell interactions are mathematically characterized, we can derive two network metrics to describe the interaction dynamics:

Clustering coefficient measures the degree to which network nodes are clustered together. In this paper, we consider a global clustering coefficient which is based on triadic network nodes 65 . Since our network links are constrained by the influential distance, a high clustering coefficient means that the network nodes are not only highly active but also in close proximity to one another, without being separated by non-active cells or other particles in the extracellular space.

Communities are defined as groups of nodes with high clustering coefficient (i.e., densely interconnected) that are only sparsely connected with the rest of the network. Therefore, a cell group may have only a few links across another group, but if the two cell groups have a large impact on one another (e.g, gene expression synchronization.), then they are considered to belong to the same community.

We also define three quantitative parameters:

QS-regulated expression: As indicated in our QS activity based network definition, the percentage of networked nodes (i.e. QS-activated cells) over the total number of cells (i.e., the sum of QS + cells, QS − cells, and QSI -resistant cells) can be used as a measure of the overall activity of QS-regulated product expression in a biofilm. Given that QS systems commonly up-regulate virulence factors, its expression is generally considered to be positively correlated with the percentage of networked nodes.

QS + genotype ratio: The ratio of the sum of the number of QS + cells and QSI -resistant cells over the total number of cells. This measure is critical because it indicates how strong the QS activity in the bacteria populations would be once the QSI treatment is reduced or stopped, regardless of how much it is inhibited when QSI is used.

Resistance: The ratio of the number of QSI-resistant cells over the total number of cells (i.e., the sum of QS+ cells, QS− cells, and QSI-resistant cells). Note that QSI resistance is defined as the continuance of QS activity following QSI treatment.

Simulation Environment Configuration

All the scenarios we model bacterial growth in a 3D microfluidic environment (250 μm × 205 μm × 400 μm) is initialized and inoculated with 25 wild-type QS+ and 25 QS− mutant cells, all of which are non-overlapping and randomly attached to the substrate, see Fig. 1. A constant nutrition concentration, S = 1 mM, available to both types of cells is maintained in the bulk layer and free to diffuse to the biomass layer. The duration of the growth period is 5 days.