Numerical investigation of microbial quorum sensing under various flow conditions

Introduction

Microorganisms preferentially reside on solid surfaces, which often leads to a closer proximity of neighboring cells than when in a planktonic form (Donné & Dewilde, 2015). At elevated cell densities, microorganisms need to efficiently coordinate the expression of energetically expensive phenotypes, such as biofilm development, exoenzyme production and microbial dispersal. Efficiency is achieved by producing and detecting relatively cheap signaling molecules which regulate the phenotype expression only when a sufficient signal concentration has been reached (Miller & Bassler, 2001). This microbial decision-making process called “quorum sensing (QS)” was originally understood as a cell-to-cell communication to identify conspecific population density and accomplish cooperative behaviors (Fuqua, Winans & Greenberg, 1994). However, a number of studies have indicated that QS is not necessarily a social trait (Redfield, 2002; Carnes et al., 2010) and depends not only on the population but also on the spatial distribution of microbial cells (Alberghini et al., 2009; Gao et al., 2016). These observations led to an alternative QS concept in which QS depends strictly on the local concentration of signaling molecules (Hense et al., 2007; Hense & Schuster, 2015). This suggests that, to understand QS processes, an integrative approach is required analyzing a multitude of factors including microbial density (Fuqua, Winans & Greenberg, 1994), production and decay kinetics (Lee et al., 2002; Fekete et al., 2010), and transport of signaling molecules through advection and diffusion (Redfield, 2002), as well as the spatial distribution of microorganisms (Alberghini et al., 2009). Thus, spatial constraints and responses may be as important as other biological considerations for the evolution and maintenance of QS. This idea is known to be true in biofilms where cooperative strategies are able to evolve if cooperators are spatially aggregated (Xavier & Foster, 2007).

Individual microbial cells synthesize and release signaling molecules at a basal rate. At low population densities, the concentration of signaling molecules remains low as it degrades both biotically and abiotically (Lee et al., 2002; Yates et al., 2002). At a sufficiently high microbial population density, however, the extracellular concentration of signaling molecules reaches a threshold concentration that activates gene and phenotypes expression (Hense & Schuster, 2015). When QS regulates the production of costly public goods, this balances production cost and the overall benefit (Pai, Tanouchi & You, 2012; Heilmann, Krishna & Kerr, 2015; Schluter et al., 2016), while under nutrient limited conditions, QS can regulate microbial dispersal (Solano, Echeverz & Lasa, 2014; Boyle et al., 2015), improving chances of survival. QS induction also often upregulates genes controlling production of signaling molecules resulting in enhanced signal production (Ward et al., 2001; Fekete et al., 2010; Pérez-Velázquez et al., 2015). Such autoinduction has been thought to confer evolutionary stability and fitness advantages (Brandman et al., 2005; Mitrophanov, Hadley & Groisman, 2010; Gao & Stock, 2018), but its effects on neighboring microbial aggregates and evolutionary benefits in a spatial context have not been fully understood.

QS induction is affected by mass transport characteristics controlling the spatial distribution of signaling molecules. In a confined space, even a single microbial cell can be QS induced if the signaling molecules accumulate to sufficiently high concentration (Carnes et al., 2010). However, higher population densities are required for QS induction in a large open space because the signaling molecules are diluted due to diffusive loss to the surrounding medium (Alberghini et al., 2009; Trovato et al., 2014). Advection may dilute the signaling molecules more effectively than diffusion and repress QS induction. Experimental observations have shown that fast advective flows increase the amount of biomass required for QS induction (Kirisits et al., 2007) and repress QS dependent gene expression (Meyer et al., 2012). Under slower flow conditions, bacteria trapped in a 3D permeable flow cell show more QS dependent gene expressions (Connell et al., 2010). QS induction can be promoted if strong advection is decoupled by heterogeneous pore geometry (e.g., dead-end pores), allowing signaling molecules to accumulate (Kim et al., 2016; Ribbe & Maier, 2016).

The signaling molecules transported either via advection or diffusion can induce QS in neighboring cells (Frederick et al., 2010; Pérez-Velázquez, Gölgeli & García-Contreras, 2016). Because the signal concentration decreases with distance from its source, cells should be located close to each other in order to send and receive enough signaling molecules to and from their neighbors (Hense et al., 2007; Matur et al., 2015). The distance between two QS induced microbial cells or aggregates is referred to as the “calling distance” and has been reported to be 5–78 μm between individual cells (Gantner et al., 2006) and ~180 μm between microbial aggregates (Darch et al., 2018). However, the dependance of QS processes on advection and diffusion suggests that transport regimes affect calling distances, highlighting the importance of relative positioning of microorganisms coupled with the mass transport characteristics of a habitat.

Here, we evaluate the effect of combined diffusive and advective transport on QS processes in environmentally relevant conditions using a reactive transport modeling approach. The advection–diffusion–reaction equation was nondimensionalized to capture the characteristic properties of QS systems (i.e., production rates of signaling molecules, cell density, mass transport and spatial distribution of microbial aggregates) and used to formulate empirical expressions describing concentration profiles of signaling molecules under various flow conditions. Using these relationships, we evaluate calling distances and threshold biochemical conditions for QS induction of a single microbial aggregate under various flow conditions. Then, we investigate QS interactions between heterogeneously distributed microbial aggregates. Finally, we demonstrate the importance of autoinduction for coordinated microbial behaviors in advection-dominated environments. This study quantifies the effect of flow velocities, autoinduction, and relative position of microbial aggregates to calling distances in a 2D flow channel.

Materials and Methods

We used the Lattice Boltzmann (LB) method to implement a numerical model for the transport of signaling molecules due to diffusion and advection. The LB method is a mesoscopic approach solving the Boltzmann equation across a defined set of particles which recovers the macroscopic Navier–Stokes equation (NSE) and advection–diffusion–reaction equation (ADRE) (Tang et al., 2013; Krüger et al., 2017). First, we obtained the flow field by solving the particle distribution function f: (1) $$f_i(\mathbf{r}+c_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)=f_i(\mathbf{r},t)+\frac{\mathrm{\Delta }t}{\tau }(f_i^{eq}(\mathbf{r},t)-f_i(\mathbf{r},t))$$where particles f_i(r, t) travel in the direction i with the lattice velocity c_i (c₀ = (0, 0), c₁ = (1, 0), c₂ = (0, 1), c₃ = (−1, 0), c₄ = (0, −1), c₅ = (1, 1), c₆ = (−1, 1), c₇ = (−1, −1), c₈ = (1, −1)) to a new position r + c_iΔt after a time step Δt. The relaxation time (τ) was described by the commonly used Bhatnagar–Gross–Krook collision operator (Bhatnagar, Gross & Krook, 1954) and the D2Q9 lattice with the corresponding equilibrium distribution function: (2) $$f_i^{eq}(\mathbf{r},t)={\omega }_i\rho (1+\frac{\mathbf{u}\cdot {\mathbf{c}}_i}{c_s^2}+\frac{(\mathbf{u}\cdot {\mathbf{c}}_i^2)}{2c_s^4}-\frac{\mathbf{u}\cdot \mathbf{u}}{2c_s^2})$$where ω_i are lattice weights (ω₀ = 4/9, ω_1–4 = 1/9, ω_5–8 = 1/36), c_s is a lattice dependent constant (here, $c_s^2$ = 1/3), and u is the macroscopic flow velocity. The moments of the discretized mesoscopic particles retrieve the macroscopic density $\rho =\sum f_i$ and momentum $\rho \mathbf{u}=\sum c_i{f}_i$. The Chapman-Enskog expansion showed that this LB approach recovers the incompressible NSE with the viscosity $\nu =c_s^2\left(\tau -{\displaystyle \frac{\mathrm{\Delta }t}{2}}\right)$ (Krüger et al., 2017). Once the flow field was obtained, we simulated solute transport with a particle distribution function g, using the regularized LB algorithm (RLB) for numerical accuracy (Latt & Chopard, 2006; Latt, 2007) and the D2Q5 lattice for numerical efficiency (Li, Mei & Klausner, 2017): (3) $${\mathrm{g}}_i\left(\mathbf{r}+{\mathbf{c}}_i\mathrm{\Delta }t,t+\mathrm{\Delta }t\right)={\mathrm{g}}_i^{eq}\left(\mathbf{r},t\right)+\left(1-{\displaystyle \frac{\mathrm{\Delta }t}{\tau }}\right){\displaystyle \frac{{\omega }_i}{2c_s^4}{\mathbf{Q}}_{\mathbf{i}}:{\mathbf{\Pi }}_i^{neq}+{\mathrm{\Omega }}_i^{\mathrm{R}\mathrm{X}\mathrm{N}}\left(\mathbf{r},t\right)}$$where c_i are the lattice velocities (c₀ = (0, 0), c₁ = (1, 0), c₂ = (0, 1), c₃ = (−1, 0), c₄ = (0, −1), c₅ = (1,1)) corresponding to the lattice weights ω_i (ω₀ = 1/3, ω_1–4 = 1/6), and Q_i : ${\mathbf{\Pi }}_i^{neq}$ is the tensor contraction of the two tensors ${\mathbf{Q}}_i={\mathit{c}}_i\cdot {\mathit{c}}_i^T-c_s^2\mathbf{I}$ and ${\mathbf{\Pi }}_i^{neq}=\sum _j{\mathit{c}}_{\mathit{i}}\cdot {\mathit{c}}_{\mathit{i}}^{\mathit{T}}\left({\mathbf{g}}_{\mathit{j}}\left(\mathbf{r}\mathbf{,}\mathbf{t}\right)\mathbf{-}{\mathbf{g}}_{\mathit{j}}^{\mathit{e}\mathit{q}}\left(\mathbf{r}\mathbf{,}\mathit{t}\right)\right)$. The reaction term in the Eq. (3) describes the production of signaling molecules: (4) $${\mathrm{\Omega }}_i^{\mathrm{R}\mathrm{X}\mathrm{N}}\left(\mathbf{r},t\right)=\mathrm{\Delta }t{\omega }_i\left(1+FH\left[\hat{A}-\hat{\theta }\right]\right)\hat{k}\hat{B}$$where F represents a multiplication factor which was set to either 0 or 10 to reflect the magnitude of autoinduced signal production (Fekete et al., 2010), $\hat{A}$ is a concentration of signaling molecules, $\hat{\theta }$ is the QS induction threshold, $\hat{k}$ is the basal production rate constant of signaling molecules, and $\hat{B}$ is the microbial density. QS induction often displays a switch-like behavior (Fujimoto & Sawai, 2013; Heilmann, Krishna & Kerr, 2015; Hense & Schuster, 2015), which is represented in the model by a step function with a higher signal production rate above the threshold concentration of signaling molecule: (5) $$H\left[\hat{A}-\hat{\theta }\right]=\{\begin{array}{c}1,\left(\hat{A}\ge \hat{\theta }\right)\\ 0,\left(\hat{A}<\hat{\theta }\right)\end{array}$$

With the imposed flow field from Eq. (1), the LB transport solver (Eq. 3) recovers the following ADRE: (6) $${\displaystyle \frac{\mathrm{\partial }\hat{A}}{\mathrm{\partial }\hat{t}}+\hat{\mathbf{u}}\cdot \hat{\mathrm{\nabla }}\hat{A}=\hat{D}{\hat{\mathrm{\nabla }}}^2\hat{A}+\left(1+FH\left[\hat{A}-\hat{\theta }\right]\right)\hat{k}\hat{B}}$$with the molecular diffusivity $\hat{D}=c_s^2\left(\tau -{\displaystyle \frac{\mathrm{\Delta }t}{2}}\right)$. Note that we are ignoring the breakdown of signaling molecules (Lee et al., 2002), limiting us to settings where production and transport are the dominant processes.

To describe the characteristic properties of a microbial system across various flow and reaction conditions, Eq. (6) was recast by introducing dimensionless quantities $A={\displaystyle \frac{\hat{A}}{\hat{\theta }}}$, $t={\displaystyle \frac{\hat{D}\hat{t}}{{\hat{l}}^2}}$, $\mathrm{\nabla }=\hat{\mathrm{\nabla }}\hat{l}$, $B={\displaystyle \frac{\hat{B}}{{\hat{B}}_{\theta }}}$, $\mathbf{u}={\displaystyle \frac{\hat{\mathbf{u}}}{\hat{U}}}$, where $\hat{l}$ is a characteristic length (i.e., the width of the flow channel), $\hat{U}$ is a characteristic fluid velocity (here, the average pore fluid velocity), and ${\hat{B}}_{\theta }$ is a threshold biomass density required for QS induction, resulting in: (7) $${\displaystyle \frac{\mathrm{\partial }A}{\mathrm{\partial }t}+Pe\mathbf{u}\cdot \mathrm{\nabla }A={\mathrm{\nabla }}^{2}A+Da}$$

This nondimensionalized ADRE is fully characterized by the Péclet number, expressing the magnitude of advective flow relative to diffusion $(Pe={\displaystyle \frac{\hat{U}\hat{l}}{\hat{D}})}$, and the diffusive Damköhler number, comparing reaction to diffusion $(Da={\displaystyle \frac{k^{\mathrm{\prime }}B{\hat{l}}^2}{\hat{\theta }\hat{D}}}$; where $k^{\mathrm{\prime }}=\left(1+FH\left[\hat{A}-\hat{\theta }\right]\right)\hat{k}{\hat{B}}_{\theta })$. A system with high Da, either due to high k’ (i.e., fast signal production), high B (i.e., high microbial density), or low $\hat{\theta }$ (i.e., high sensitivity to signaling molecules) - is more likely to be QS induced.

An important property of Eq. (7) is that its solution linearly scales in Da (Lin, Xu & Zhang, 2020). For example, if Da is increased 2-fold at a fixed Pe condition, the concentrations of signaling molecule are doubled. This linearity allows to calculate the concentration distribution of signaling molecules for any Da from a single simulation result with an arbitrary Da at a given Pe. However, this simple approach cannot be applied to the flow conditions because the solution is not linear in Pe. Therefore, multiple numerical simulations were carried out with 24 Pe conditions (Pe ∈ {0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10}) while Da was fixed at 5. For the 2D simulations in a straight channel (Fig. 1), the flow field was established by imposing pressures at in- and outlet and no flow conditions at the top and bottom boundaries, resulting in a flow from left to right. Fixed concentration (A|_{left boundary, x=0} = 0) and no-gradient $(\mathrm{\partial }A/\mathrm{\partial }x{|}_{\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{a}\mathrm{r}\mathrm{y},\mathrm{x}=4}=0)$ boundary conditions were imposed at the inlet and outlet boundaries, with no-flux at the top and bottom boundaries, respectively. All simulations were run to steady state.

Simulations were conducted for a 2D flow channel of non-dimensional length of 4 and a width of 2, discretized with 2,000 × 1,000 grid elements. The flow field (Eq. 1) was generated by imposing fixed pressures at inlet (x = 0) and outlet (x = 4) with no flow boundaries in both normal and tangential direction at the bottom (y = 0) and top (y = 2) of the domain resulting in parabolic Poiseuille flows. Simulations were carried out under low Mach numbers (Ma = u/c_s ≪ 1) to ensure incompressible flow conditions (Krüger et al., 2017).

Results and Discussion

QS processes of a single microbial aggregate

The effect of various flow conditions on the distribution of signaling molecules (A) produced from a single microbial aggregate assumed a source constrained to a single grid cell located at x = 1 was investigated under various Pe conditions (0.5 ≤ Pe ≤ 10) while Da was fixed at 5 (Fig. 1). The environmentally relevant range of Pe was chosen (Battiato et al., 2009, 2011) while Da is arbitrary because of the linearity of Eq. (7) in Da. The QS induction enhancing the signal production rate was not considered.

The signal concentration fields developed under various advective flows show maximum concentrations (A_max = A(x=1)) decreasing with increasing Pe (i.e., faster advective flow): A_max decreased from 1.68 (Pe = 1) to 1.35 (Pe = 5) and 1.21 (Pe = 10). However, A_max of all of the simulations with Da = 5 exceeded 1 (i.e., $\hat{A}$ ≥ $\hat{\theta }$), indicating the potential for QS induction. The threshold Da for QS induction (Da_θ), where A_max = 1, can easily be computed using the linearity of the nondimensionalized ADRE in Da (Eq. 7). For example, Da_θ at Pe = 1 was calculated by dividing Da = 5 by its corresponding A_max = 1.68 which resulted in Da_θ = 2.98. Thus, at Pe = 1, conditions for which Da ≥ 2.98 lead to or exceed the concentration of signaling molecules needed for QS induction. Figure 2 shows the calculated Da_θ for each simulated Pe condition.

The regression analysis revealed that the simulated Da_θ for QS induction varies as a function of Pe following the power law: (8) $$D{a}_{\theta }=1.3812P{e}^{0.2626}+1.592$$

The increasing Da_θ along with the increasing Pe indicates higher Da (i.e., higher microbial density (B), higher signal production rate constant (k’), or lower QS induction threshold ($\hat{\theta }$)) is required for QS induction under higher Pe. This result corresponds to the observed repressed QS induction under the presence of advection (Vaughan, Smith & Chopp, 2010; Meyer et al., 2012; Kim et al., 2016) and matches the pattern of biomass required for QS under varying flow conditions (Kirisits et al., 2007). Equation (8) was further evaluated by applying the experimentally measured QS parameters of Pseudomonas putida ($\hat{k}$ = 2.3 × 10⁻¹⁰ nmol/cell/h, and $\hat{\theta }$ = 70 nmol/L (Fekete et al., 2010)) in a flow system where $\hat{l}$ = 1 cm and $\hat{D}$ = 3.0 × 10⁻¹⁰ m²/s (Dilanji et al., 2012). Our results show ${\hat{B}}_{\theta }$ of 9.77, 12.2 and 13.5 × 10⁶ cells/mL at Pe = 1, 5 and 10, respectively. If Eq. (8) is extrapolated to diffusion only transport condition (Pe = 0, Da_θ = 1.592), ${\hat{B}}_{\theta }$ is estimated as 5.23 × 10⁶ cells/mL which largely agree with the experimental observation of 2.69~6.23 × 10⁶ cells/mL where signal concentration starts to show a strong spike (Table S1 in Fekete et al., 2010).

In addition to reducing A_max, advection also influenced the spatial distribution of signaling molecules. We define the “transport distance” (d) as the distance between the point of production (x₀) and the point (x₁) where the signal concentration reaches a certain value A* (i.e., d = |x₀ − x₁|), distinguishing it from the “calling distance” between two QS induced microbial cells or aggregates. If the signal transport occurred only through diffusion, transport distances would be isotropic (Alberghini et al., 2009). However, advection resulted in anisotropic concentration distribution where upstream transport distances (d_up) are much shorter than the downstream distances (d_dn). Moreover, fast advective flows (i.e., high Pe) reduced overall transport distances which are illustrated in Figs. 1A–1C as the shrinking areas covered by contour lines. For example, the (nondimensional) transport distances to the location where A = 0.1 are d_up = 0.08 and d_dn = 0.62 at Pe = 1 and Da = 5 (Fig. 1A). These values decrease to d_up = 0.033 and d_dn = 0.27 at Pe = 5 (Fig. 1B) and to d_up = 0.023 and d_dn = 0.19 at Pe = 10 (Fig. 1C).

Empirical approximation of concentration profiles

Obtaining transport distances for different Pe conditions requires running numerical simulations for each of the corresponding Pe. However, this may be avoided if we can express the concentration profiles as a function of Pe. For this purpose, parametric regression analysis was applied to the numerically obtained concentration profiles along the bottom of the flow channel (Fig. 3).

Several parametric regression models (linear, power, exponential and polynomial models) were tested to the upstream (A_up(x); 0 ≤ x ≤ 1) and downstream (A_dn(x); 1 < x ≤ 4) signal concentration profiles. Among the tested regression models, the exponential (Eq. 9) and power-law models (Eq. 10) provided the best fit for log-transformed upstream and downstream signal concentration profiles, respectively. In the regression analysis of upstream profiles, only the locations where A(x) > 0.001 were used to improve the fitting quality and the signal concentration at x = 1 was fixed as 1. The additional regression analysis was then carried out for the coefficients (a, b and c) obtained from simulated profiles at 24 Pe conditions to construct a relationship between the coefficients and Pe (Figs. 3B–3D). The exponential and power-law models provided the best fit for the upstream (Eqs. 11–13) and downstream coefficients (Eqs. 14–16), respectively: (9) $$A_{\mathrm{u}\mathrm{p}}\left(x\right){|}_{x\le 1}=\exp \left(a_{\mathrm{u}\mathrm{p}}\left(x^{b_{\mathrm{u}\mathrm{p}}}-x^{c_{\mathrm{u}\mathrm{p}}}\right)\right)$$ (10) $$A_{\mathrm{d}\mathrm{n}}\left(x\right){|}_{x>1}=\exp \left(a_{\mathrm{d}\mathrm{n}}\ln {\left(x\right)}^{b_{\mathrm{d}\mathrm{n}}}+c_{\mathrm{d}\mathrm{n}}\right)$$where A_up and A_dn are 0 in the down- and up-stream directions, respectively, and (11) $$a_{\mathrm{u}\mathrm{p}}=0.376\exp \left(-2.5975Pe\right)+2.7165\exp \left(0.0244Pe\right)$$

(12) $$b_{\mathrm{u}\mathrm{p}}=20.311\exp \left(0.0115Pe\right)-13.38\exp \left(-0.6121Pe\right)$$ (13) $$c_{\mathrm{u}\mathrm{p}}=-7.1289\exp \left(0.0348Pe\right)+5.9469\exp \left(-0.4272Pe\right)$$ (14) $${𝑎}_{\mathrm{d}\mathrm{n}}=8.6156𝑃𝑒-0.0668-13.3056$$ (15) $$b_{\mathrm{d}\mathrm{n}}=0.1051P{e}^{-0.2522}+0.1082$$ (16) $$c_{\mathrm{d}\mathrm{n}}=-7.5322P{e}^{-0.0464}+8.7195$$

Equations (9) and (10) can be used as approximations of the concentration profiles along a pore channel without running simulations for various Pe conditions, with the microbial aggregate located at x = 1. Due to the linearity in Da, the concentration profiles at different Da conditions can be calculated simply by multiplying Da/Da_θ to Eqs. (9) and (10), so that (17) $$A\left(x\right)={\displaystyle \frac{Da}{D{a}_{\theta }}\left(A_{\mathrm{u}\mathrm{p}}\left(x\right)+A_{\mathrm{d}\mathrm{n}}\left(x\right)\right)}$$

These analytical expressions are applicable not only to QS but also to other chemical processes subject to zero-order production reactions (e.g., Bezemer et al., 2000; Tang et al., 2015). The equations become less accurate at low Pe as under low flow conditions, the estimates from Eq. (17) in a flow channel with a small width (i.e., low $\hat{l}$ and Pe) could underestimate the actual concentration because the confined channel width would push the produced chemical further upstream and downstream.

The effect of QS induced signal production on transport distances

QS often involves autoinduction which substantially increases signal production rates. The effect of autoinduction on transport distances was investigated by using Eq. (17) for the conditions without (F = 0; Da = Da_θ) and with (F = 10; Da = 11Da_θ) enhanced signal production. The transport distances from a single microbial aggregate under various Pe were then calculated using Eq. (17) for the location x.

Figure 4 shows the transport distances without (Fig. 4A) and with (Figs. 4B and 4C) the enhanced signal production at Pe = 1, 5 and 10. The concentration ratios (0.1 ≤ A/A_max ≤ 0.9) were used instead of absolute concentrations to generalize transport distances for various Da conditions. For example, the transport distance (d_A) for A/A_max = 0.5 indicates that A(x₀ + d_A) = 0.5 if Da = Da_θ while A(x₀ + d_A) = 0.05 when Da = 0.1Da_θ. The consequence of the enhanced signal production was the significant increase of d_up and d_dn. Without the enhanced signal production, d_up and d_dn for A/A_max = 0.4 at Pe = 1 were estimated as 0.021 and 0.024, respectively (Fig. 4A). These values increased to d_up = 0.1 and d_dn = 1.28 with the enhanced signal production (Figs. 4B and 4C). The downstream transport distance of 1.28 is translated into 6.4 mm in a flow channel with $\hat{l}$ = 1 cm. This result is much longer than the generally observed ranges of calling distances (Whiteley, Diggle & Greenberg, 2017). However, we emphasize again that the transport distance merely indicates the distance of signaling molecules transported from a source location while the calling distance involves QS induced microbial cells or aggregates.

QS induction between spatially distributed multiple microbial aggregates

QS processes of multiple aggregates were investigated by constructing the concentration profiles using Eq. (18). Concentration fields of signaling molecules with multiple microbial aggregates can be calculated as the superposition of the concentration profile produced by each individual aggregate: (18) $$A\left(x\right)=\sum _{i=1}^n{\displaystyle \frac{D{a}_i}{D{a}_{\theta }}\left(A_{\mathrm{u}\mathrm{p}}\left(x+d_{i0}\right)+A_{\mathrm{d}\mathrm{n}}\left(x+d_{i0}\right)\right)}$$where n is the number of aggregates, d_i0 is the distance between x_i and x₀ (d_i0 = x_i − x₀), x_i is the location of ith aggregate, x₀ is the reference location (x₀ = 1), Da_i is the Da calculated only with the density of ith microbial aggregate (i.e., microscopic Da), and A_up and A_dn are Eqs. (9) and (10), respectively. Here, an example system with macroscopic $Da(D{a}_T=\sum D{a}_i)$ = 3.2Da_θ consist of four aggregates ($\mathbb{A}$_1–4) located at x₁ = 0.4, x₂ = 1, x₃ = 1.096 and x₄ = 1.7 with the evenly distributed microscopic Da_i (i.e. Da₁ = Da₂ = Da₃ = Da₄ = 0.8Da_θ) was tested. In using Eq. (18), the profile was first constructed for Da_i = Da* that does not consider autoinduction (F = 0). Then, if there is an aggregate with A(x_i) ≥ 1, the profile was recalculated with updated Da_i = (1 + F) × Da* until all Da_i with A ≥ 1 were updated.

The signal concentration profile produced by four aggregates without the enhanced signal production (F = 0) illustrates the crucial importance of relative positioning of microbial aggregates for QS induction with respect not only to each aggregate but also to the flow direction (Fig. 5A). The microscopic Da_i was set such that the maximum concentration produced by a single aggregate was 0.8, as observed at the most upstream location ($\mathbb{A}$₁ at x₁ = 0.4). But due to transport, the local concentration at $\mathbb{A}$₂ reached 0.879, receiving A of 0.048 and 0.031 from $\mathbb{A}$₁ and $\mathbb{A}$₃, respectively. $\mathbb{A}$₃ received slightly less signaling molecules from $\mathbb{A}$₁ (A = 0.044) due to the longer distance of $\mathbb{A}$₃ than $\mathbb{A}$₂ from $\mathbb{A}$₁. However, $\mathbb{A}$₂ provided much more signaling molecules (A = 0.157) to $\mathbb{A}$₃ than was provided by $\mathbb{A}$₃ because of advective flows favoring downstream transport of signaling molecules (Figs. 2 and 4). As a consequence, $\mathbb{A}$₃ exceeded the QS threshold (A(x₃) = 0.044 from $\mathbb{A}$₁ + 0.157 from $\mathbb{A}$₂ + 0.8 from $\mathbb{A}$₃ + 0 from $\mathbb{A}$₄ = 1.001 > 1) while the upstream located $\mathbb{A}$₂ did not. The QS induction of $\mathbb{A}$₃ demonstrates the importance of transport distances. QS induction was achieved because of the upstream aggregates located within the transport distance of 0.696. However, the calling distance would have been estimated as the length of a grid voxel (0.002) because only $\mathbb{A}$₃ was QS induced. Therefore, considering only the calling distance could lead to the wrong conclusion that the local Da condition at $\mathbb{A}$₃ (i.e., Da₃ = 0.8Da_θ) is a sufficient condition for QS induction. Although $\mathbb{A}$₄ did not reach the QS induction threshold, it received A from all the other aggregates resulting in a concentration (A(x₄) = 0.029 + 0.044 + 0.048 + 0.8 = 0.921) that was higher than at $\mathbb{A}$₂ despite the longest separation distance from other aggregates.

Accounting for QS induction (F = 10) increased the transport distances and hence induced other aggregates (Fig. 5b). With the same spatial distribution, QS-induced $\mathbb{A}$₃ produced signaling molecules much more and faster (i.e. k’ = 11$\hat{k}$ and Da₃ = 8.8Da_θ) and provided more signaling molecules to $\mathbb{A}$₂. As a result, A(x₂) exceeded the QS threshold (0.048 + 0.8 + 0.335 + 0 = 1.183). The QS induction of $\mathbb{A}$₂ and $\mathbb{A}$₃ resulted in the final signal concentrations of A(x₂) = 9.183 (= 0.048 + 8.8 + 0.035 + 0) and A(x₃) = 10.569 (= 0.044 + 1.725 + 8.8 + 0). While $\mathbb{A}$₄ still did not contribute signaling molecules to any of upstream aggregates, enhanced contribution from $\mathbb{A}$₂ and $\mathbb{A}$₃ QS induced $\mathbb{A}$₄, A(x₄) = 9.839 (= 0.029 + 0.48 + 0.53 + 8.8). Despite increased transport distances by QS induction, $\mathbb{A}$₁ was still too far away from the other aggregates thus the signal concentration at $\mathbb{A}$₁ remained unchanged A(x₁) = 0.8. As a result of the QS induction of $\mathbb{A}$_2-4, Da_T had increased from the initial 3.2Da_θ (= 0.8Da_θ × 4) to 27.2Da_θ (= 0.8Da_θ + 3 × 11 × 0.8Da_θ).

This example illustrates the importance of enhanced signal production on the spatial propagation of QS induction. While only $\mathbb{A}$₃ experienced signaling molecule levels that could induce QS when all the aggregates produce signaling molecules at the basal production rate, the enhanced signal production of $\mathbb{A}$₃ when considering induced production (F = 10) provided more signaling molecules to its adjacent microbial aggregates and resulted in the QS induction of neighboring aggregates, $\mathbb{A}$₂ and $\mathbb{A}$₄. It may be counterintuitive that the upstream-located $\mathbb{A}$₂ was also QS-induced by the contribution from $\mathbb{A}$₃ despite the contracted upstream transport distances under the presence of advective flows. This result shows that the enhanced signal production can overcome the influence of advection and promote QS induction, and provide a way to provoke upstream microbial aggregates, for example, to slow down the substrate consumption to ensure efficient resource utilization in crowded environments (An et al., 2014).

Conclusions and Perspectives

This study has demonstrated that advection and the enhanced signal production can determine the spatial extent of QS induction. Reactive transport simulation results reveal that fast flow conditions dilute signaling molecules and thus higher Da_θ (i.e., faster signal production rate, higher microbial density, or lower QS induction threshold concentration) is required for QS induction. Reduced upstream delivery of signaling molecules under advective flow limits propagation of QS; it can be relaxed if autoinduction increases signal production rates. Our study results highlight the importance of relative positioning of microbial aggregates with respect to flow directions and the role of autoinduction to overcome advection for upstream signal propagation.

The simulations focused on the effect of various flow conditions on QS and assumed that microbial aggregates have a negligible impact on flow fields, which is a reasonable approximation for low microbial density conditions. However, it may not hold when large aggregates producing extracellular polymeric substances (EPS) perturb flows substantially. In such a case, estimating signal transport requires fully resolving nonlinear feedback between cell activity and fluid flow (Thullner & Baveye, 2008; Carrel et al., 2018; Jung & Meile, 2019), including diffusion limitation (Stewart, 2003). Finally, accounting for degradation of signaling molecules (Lee et al., 2002; Yates et al., 2002) and increased spreading of signaling molecules in 3D systems than 2D, would result in shorter transport distances than this study.

Although QS mediated gene expression has been understood as evolutionarily beneficial collective behaviors, long transport distances observed in this study suggests that it may not be always true. The transport of signaling molecules, especially in downstream direction, combined with enhanced signal production, suggests that QS induction can be decoupled from microbial density. In the above example (Fig. 5B), any microbial cell located where A > 1 (e.g., A(x = 2.5) = 1.05) would have been QS-induced, independent of the local cell density. This could lead to detrimental impacts on a microbial population, unless there are other counteracting mechanisms such as differential QS induction sensitivity to signal concentration even in within a clonal population (Darch et al., 2018) or biofilm formation modifying local transport characteristics (Emge et al., 2016). Future investigations should explicitly examine the evolutionary consequences of QS strategies in spatially heterogeneous environments under advective–diffusion–reaction dynamics.