Burst statistics in an early biofilm quorum sensing model: the role of spatial colony-growth heterogeneity

Quorum-sensing bacteria in a growing colony of cells send out signalling molecules (so-called “autoinducers”) and themselves sense the autoinducer concentration in their vicinity. Once—due to increased local cell density inside a “cluster” of the growing colony—the concentration of autoinducers exceeds a threshold value, cells in this clusters get “induced” into a communal, multi-cell biofilm-forming mode in a cluster-wide burst event. We analyse quantitatively the influence of spatial disorder, the local heterogeneity of the spatial distribution of cells in the colony, and additional physical parameters such as the autoinducer signal range on the induction dynamics of the cell colony. Spatial inhomogeneity with higher local cell concentrations in clusters leads to earlier but more localised induction events, while homogeneous distributions lead to comparatively delayed but more concerted induction of the cell colony, and, thus, a behaviour close to the mean-field dynamics. We quantify the induction dynamics with quantifiers such as the time series of induction events and burst sizes, the grouping into induction families, and the mean autoinducer concentration levels. Consequences for different scenarios of biofilm growth are discussed, providing possible cues for biofilm control in both health care and biotechnology.


System and Numerical Methods
Certain bacteria are known to detect foreign microorganisms or interpret multiple signals via the mechanism of QS [49][50][51] . We here focus on the mechanism of detecting their own species with a single signalling molecule species called the AI. This name stems from the fact that the emitting cell senses its own signalling molecule, and the extracellular AI concentration feeds back on the cell's gene-regulation cascades 14,48 .
We note that for such intra-species communication Gram-negative bacteria such as Vibrio fischeri or Pseudomonas aeruginosa produce species-specific molecules from the class of N-acyl homoserine lactones (AHL) as transmitters 20,52,53 . In total 56 different AHLs were identified so far 54 . The bacterial cells use specialised oligopeptide receptor types to measure the extracellular concentration of AHL molecules 15 . After secretion the AI molecules act as a passively diffusing signal spreading until they either degrade or bind to and are detected by a receptor. The AHL-driven regulatory QS circuits are mainly governed by the Luxl-type protein responsible for AHL synthase, and the Luxr-type protein involved in AHL-response regulation 52 . These proteins act as transcription factors regulating the expression of QS-related genes. Below a given extracellular concentration of AHL molecules the AI synthase is expressed at a lower, basal level. However, above a certain threshold, in a positive feedback loop the cell changes the expression level of the gene controlling the synthase of the AHL molecules. This process is called "induction" and is a direct response to the local cell concentration.
Since every QS-responsive cell in a single-species environment is secreting and detecting these molecules, the concentration-increase effect can build up quickly when critical population densities are reached. Biofilm growth typically starts on a surface, only at later stages biofilms grow into the third dimension. The classical experimental scenario is a biofilm grown in a Petri dish, with the vertical height of the solution much smaller than the horizontal extension (see also Fig. 1). We choose below a two-dimensional description; the formulation can be extended to the third dimension (on top of the growing colony). We use a protocol of simulations written in R-code. The physical parameters classifying our system is the typical cell radius r c and the initial number n c (t = 0) of seed cells in the circular simulations area of radius R.
The cell division rate is γ, and-as a simple model of bacterial motility or diffusion-we place the two daughter cells such that one next-generation cell occupies the same position as the mother cell {x old , y old }, while (according to the first modelling scenario) the second cell is displaced by a fixed distance and at a random azimuthal angle in the plane. In the second modelling scenario, the positions of displaced daughter cells obey a Gaussian distribution with the standard deviation σ str along x and y axes, namely, We checked that both cell-spreading scenarios yield qualitatively similar results. In the figures below we provide the values of d new and σ str for the situations when, respectively, the first or the second scenarios were implemented in the simulations. We note that this scenario is a first approximation: in reality the direction of cell placement depends on the cell population density and displacement separations should reflect bacterial mobilities. Note that a number of biofilm-forming bacteria move in reality using flagellar motion (and do not jump in space). Moreover, the emerging flows can cause dispersion and the motility-to-biofilm transition is known to exist 55 . As we here focus on the effects of geometric disorder, however, neglecting such dynamic effects appears to be in order.
AI molecules are produced with a basal production rate ϕ; after induction ϕ is enhanced by the ratio η such that induced cells produce AIs with rate ϕ(1 + η). The decay rate of AI molecules is δ and they diffuse with the diffusion coefficient D. The threshold concentration for induction of AI QS-signalling molecules is β thr .
We set the radius of the circular model cells to r c = 2.5 μm, the initial colony radius to R = 100 μm (unless specified otherwise), and the initial cell density to ρ = 0.0025 cells per μm 2 . The mean number of cells 〈n c 〉 = ρπR 2 ≈ 78.5 (or 〈n c 〉 = ρπR 2 ≈ 314 when R = 200 μm) is chosen from a Poisson distribution (see below). The cell division rate is γ = 1/[60 min], the basal AI production rate in the model is ϕ = 1/sec, and the induction threshold is β thr = 30 AI molecules per μm 2 on the surface (unless otherwise indicated). All these values are considered as fixed, biologically motivated parameters [19][20][21][22][23][24][25][26] . The variations we study reflect changes in the spatial disorder of the growing cell colony. Otherwise the parameter space would be too large for the current scope.
Different approaches in simulating the growth of the bacteria colony can be chosen, such as lattice-free or lattice-based, interacting or non-interacting, as well as individual-or density-based. In lattice-based approaches new cells can only occupy points of the lattice, while in lattice-free approaches the whole plane is available for cell placement. A detailed comparison between the two approaches is reported in 56 . We choose the more natural lattice-free strategy. We follow an individual-based approach in which we keep track of each cell, instead of continuous cell densities measured in approaches modelling the colony dynamics in terms of reaction-diffusion equations 57 . To keep the simulations sufficiently simple we also choose a non-interacting approach, such that, in principle, placed cells may overlap with already existing ones. This effect is not too pronounced, as it can be interpreted as cells growing on top of each other (the next step in volume-expanding biofilm growth) 28,58,59 .
To construct the overall concentration perceived by a single cell in the colony, we superimpose the AI fields c(r, t) produced by individual cells, as mediated by the reaction-diffusion equation 2 in polar coordinates, where r is measured from the cell centre. On top of this diffusive spreading dynamics the AI molecules degrade exponentially with rate δ, while we take the AI production rate ϕ to be unity in our time units (this choice merely renormalises time). As AI molecules degrade sufficiently fast, we choose natural boundary conditions for c(r, t) in an open system requiring the density to vanish far away, lim r→∞ c(r, t) = 0. We refer to 29 for the theoretical quantification of other-adsorbing and reflecting-boundary conditions on the properties of QS. For a colony in a confined area (e.g., a Petri dish), the AI concentration profile might get reflected from the boundaries and thereby affect the QS properties (especially for slowly degrading AIs and small confining domains). A multitude of approaches for the computational modelling of cell-cell communication and bacterial colony growth has been developed (see, e.g., references) 47,60-73 . They involve stochastic modelling on individual level, modelling approaches based on a variety of deterministic partial differential equations (they describe the transport properties of various molecules, mass fluxes of the components, diffusion and advection properties, etc.), or stochastic cellular-automata models, to mention but a few. Multiple aspects of QS were considered on different levels of sophistication, with a number of features of natural bacterial habitats being described.
www.nature.com/scientificreports www.nature.com/scientificreports/ The solution of the AI diffusion Equation (2) rather rapidly relaxes to the steady state. We thus arrive at the time-independent single-cell AI concentration field (see Fig. 2 c st 0 0 in terms of the modified Bessel function K 0 of the second kind (Macdonald function). Here we chose the integration constant such that the concentration is unity (in units of the induction-threshold level β thr ) for the midpoint distances of two cells that just touch each other. In the case of overlapping cells, when r ≤ 2r c , the concentration is truncated to the value at δ = r r D 2 / c . The diffusion coefficient D and degradation rate δ solely occur as the ratio This parameter corresponds to the inverse square of the signal range (the typical distance reached by a diffusing molecule before degradation). Finally, the overall concentration of AIs perceived by a single cell is simply the sum of the concentration fields c st (r) originating from all cells of the colony. In this superposition we assume that the steady-state concentration c st (r) is reached much faster than the cell division time 1/γ.
The half-life estimates for AHL molecules, 1/δ, vary from 10 hours 28 over 1 day 74 to 7 days 29 (in the absence of bacterial cells). However, in the case of absorbing system boundaries the depletion can be much faster: namely, of the order of the time needed for a AI molecule to diffuse across the system and get absorbed to the boundaries (some hours) 29 . These variations give rise to variations of the parameter α. Moreover, AI molecules get depleted by the QS-cells themselves and can be washed away by local fluid flows. The measured characteristic distances in the case of an absorbing boundary proposed in 28 is δ = D/ 3 mm which translates to α = δ/D ~ 0.1/mm 2 . The values of α in our model vary between 10 −3 and 10 3 /μm 2 that corresponds to the range of typical distances between ≈0.03 and ≈30 μm. The separation of 30 μm is indeed a typical extension of early-stage biofilms, both in QS experiments 14 and modelling approaches 49 , while 0.03 μm is clearly a very short length of spatial extensions for a bacterial colony. On the other hand, the boundaries in the current model are open and not absorbing. Finally, the measured diffusion constant of AHL molecules in agarose discs is D ~ 300 μm 2 /sec 28,29 . In the model we set a constant, bacteria-density-independent diffusion coefficient of AI molecules (a valid approximation for the early stages of colony development). Moreover, we refer to 75 for the absolute and relative-to-water diffusivities measured in bacterial biofilms for a number of solutes of different molecular weights. Lastly, the basal production rate of AI molecules by a single bacterial cell is in the range ~110 per sec 76 .
The initial condition is the number n c (t = 0) of seed cells, spread randomly within the fixed colony radius R. To reflect natural variations every sample colony starts with a different number of cells, that is, n c (t = 0) is a realisation of a Poisson-distributed random variable with Entering the actual division loop, the mother cell for division is picked uniformly-randomly. The spawn position of the daughter cell (the other daughter cell takes up the same position as the mother cell) is determined by picking a random angle and a fixed distance d new . Alternatively, the new position of the daughter cell can be calculated using the Gaussian distribution (1), for which the degree of dispersion is modulated by σ str . Every seed cell is assigned a label that gets inherited by all daughter cells. During the colony evolution we keep track of a cell's origin and identify "genetically related" cell families or lineages. Right: Distances r i (measured in μm from the cell centre) at which the concentration has dropped to i% of its initial value, as function of α.
The next step is to calculate the AI concentration at the new cell location, to check if this cell may be induced right away. If a cell is induced its AI synthase is up-regulated by the factor (1 + η) leading to stronger emission of AIs 47 . After the placement of the daughter cell the AI levels of all cells are updated to check if more cells are induced. If not, the increment t gap is added to the time counter, and a new mother cell is picked. This iteration continues until the designated termination condition is reached. Colonies with a high value of d new or, respectively, σ str result in somewhat larger but, importantly, more spatially homogeneous colonies due to mixing of daughter cells from different seed cells. Conversely, smaller values of these spatial homogeneity parameters lead to smaller, more densely clustered colonies with a higher degree of family separation (see below). We study the system by variation of the inverse signal-range parameter α, the spatial homogeneity parameters d new or σ str , and the surplus amount η of AI production upon induction. We expect that AI concentrations are higher in densely crowded areas and thus colonies with smaller spatial homogeneity have their first-induction point earlier. Conversely, it may take longer until the entire colony is induced because subcolonies can be isolated and thus cannot profit from high AI levels due to the patchiness when the signal range is limited. Otherwise, if η is large enough the increase in signal production of an even spatially confined induced area might, given a sufficiently large signal range, lead to a cascade inducing larger fractions of the population. To explore how far such effects indeed occur we quantify both on a single-colony level and on a statistical ensemble level.
Since measuring AI levels of every cell is unrealistic experimentally, we calculate the arithmetic mean concentration 〈 〉 c of all cells as a central object for our statistical analyses. Before studying the time evolution of 〈 〉 c for a given realisation of cell colony we consider the relation of 〈 〉 c with the signal range and the number of seed cells n c . To this end a fixed number of randomly spread seed cells within a fixed radius R are chosen. The resulting 〈 〉 c as seen in Fig. 3 then shows (i) a deterministic relation to α. (ii) The corresponding sample variance s 2 in Fig. 3 displays a similar behaviour. It describes the range of possible mean concentrations with a given signal range, that is, it helps to appraise the sensitivity of c to the exact spatial distribution of the seed cells.
While higher expected mean concentrations using configurations with smaller α values are immediately intuitive, the lower sample variance for higher α values leads to the conclusion that the concrete spatial configuration of cells is playing a minor role. Considering that large values of α imply an effective signal range of barely a few microns, as can be seen in Fig. 2, means that in this case variations are mostly caused by close-by or even overlapping cells. The extent of cell overlap and overall degree of spatial heterogeneity of a colony on the induction process is elaborated further when investigating the statistical data in section 4. Before doing so, we first evaluate the time evolution of single colonies, to see what parameters are relevant to classify the induction dynamics of a growing heterogeneous colony. Fig. 4 we show the time evolution of a growing cell colony during the time span from 90 to 240 min. We note that for different cell types or environmental conditions this time scale may need to be significantly rescaled while the effects of geometric disorder remain unaltered. In the colony of Fig. 4, due to the random placement of daughter cells denser clusters of cells emerge. At 192 min the first-induction event is recorded in the densest area, as can be seen at time 210 min (dark blue cells at the location near (−90, 40)). During the following 30 min interval several other local clusters are induced. While in this simulation the initial seed cells were placed randomly within the radius R, even for perfectly regularly placed seed cells an apparent randomisation emerges relatively quickly due to the random placement of daughter cells, see Fig. S.1. It is typical for most parameter choices used here that induction occurs after some 180 min, while a large part of the entire colony is induced after some 250 min.

Single-colony time evolution. In
As shown in Fig. 2 the AI concentration field of a single cell decays rapidly over a few μm for the chosen values of the inverse signal-range parameter α. This suggests that even around very dense areas the cumulative AI concentration has a rather short range, due to degradation of AI molecules. Figure   www.nature.com/scientificreports www.nature.com/scientificreports/ reasoning: the concentration variation depicted in the heat map for the shown bacteria colony at the moment of the first-induction event is relatively high. The concentration in areas with low cell densities is almost vanishing, in the denser regions an up to three-fold change of the concentration level can be observed. This effect allows for rather well defined, local induction regions, as opposed to an homogeneous situation as that considered in 47 . The first induction occurs after 192 min. The colony radii of ~100 μm in our simulations agree with previous estimates/evidence 6,48,49 (larger colonies are realised for higher σ str values). The cell coordinates in the panels are given in μm. Figure 5 demonstrates the effect of spatial disorder on the induction dynamics of the growing cell colony. When daughter cells are placed relatively close to their mother cells (d new = 5 μm) the initial inhomogeneities of the seed cluster are amplified, and significantly denser local clusters emerge. In these clusters induction occurs relatively quickly, see also the quantitative analyses below. Upon increasing the placement distance of the daughter cells up to d new = 25 μm, the colony expands significantly (note the different scales in the four panels) and, concurrently, appears much more homogeneous. In such a scenario it naturally takes longer for the first induction to occur. It will be an interesting question to address whether there eventually occurs an inversion towards full colony induction, that is, whether homogeneous colonies are completely induced earlier than heterogeneous ones.
The growing colony is characterised by different quantities in  www.nature.com/scientificreports www.nature.com/scientificreports/ η upon induction of a cell. How precisely this induction behaviour is influenced by the magnitude of η is investigated later on.
We first visualise the burst effect in a single-colony analysis. To see the extent of induction bursts we identify the largest jumps in c and χ corresponding to the largest peaks of k t . These are compared to the average number of induced cells per time step so far. The average is taken over a sample including all division steps since the last burst event. If the difference in burst size is large enough, the algorithm colours all cells that got induced at the critical time step. For the algorithm to activate a threshold difference of 8 standard deviations to the average jump size was chosen. Figure 6 decomposes the distribution of cell family sizes and induction levels per lineage after 240 min, further detailing the burstiness of the induction process.
Examining the sample colony from Fig. 7 it is indeed possible to visualise the induction steps. The first induction in this simulation occurs after 165 min, in a densely crowded region in the lower right quadrant near the origin (0, 0), see Fig. 7(a). The first burst (124 cells induced, red) was large enough to induce almost 10% of the entire colony at this time. Only minutes later another area with 134 cells (green) was induced close to the first region. Within less than five minutes χ almost reached 20% of the colony. Though the induction process was quickly accelerating, it slowed down right after these two bursts since the two induced clusters are spatially quite clearly separated from the rest of the colony and thus beyond signal reach. If these two regions had been more entangled they could have likely acted as a pacemaker for the induction of other larger vicinal clusters. (Note that for a region to become induced it is not only necessary to have a high local cell density but induction probabilities are increased when other induced clusters are within signal range). After 20 more minutes the next induction period heralds. In the lower left quadrant a large region comprising 140 cells (yellow) are induced. Three major bursts follow shortly after, see Fig. 7(b,c), leading to the induction of nearly the whole quadrant. This period lasts for less than 15 min and pushes χ to almost 60%. This sample colony evolution also underlines that on average induction bursts appear to decline in size with time. www.nature.com/scientificreports www.nature.com/scientificreports/ When looking at the colony after 240 min one can see that another large burst becomes more and more unlikely. After the induction hot spots are distributed over all regions they act as pacemaker for the induction of their vicinity surrounding. After 220 min of time evolution, see Fig. 7(d), and 50 min after the first induction a larger fraction of cells is coloured dark blue indicating that they were induced little by little, or in small bursts, see Fig. 7(e,f). This behaviour can be explained by the AI concentrations of many cells building up to the point of www.nature.com/scientificreports www.nature.com/scientificreports/ first induction close enough to the threshold β thr so that just a few induced cells next to or within the cluster with additional AI secretion tip the balance towards cluster induction.
A further way to analyse the induction of the colony is depicted in Fig. 8 and S.4. At different instances of the time evolution of colonies with different degree of disorder (determined by the standard deviation σ str of the Gaussian daughter cell displacements) the induction level is characterised by the histograms for different induction levels χ, as shown in Fig. 8. It can be seen that homogeneous colonies contain hardly any induced cells even at quite long times, while heterogeneous colonies already have quite advanced induction levels. However, in the course of time the homogeneous colonies quite quickly catch up, as induction occurs as a relatively sharp event, when the colony becomes critical almost everywhere. Concurrently, induction bursts increase in size with the degree of homogeneity of the colony. We show in Fig. S.4 the effect of the AI production enhancement η when a cell is induced. Larger values of η clearly speed up the induction dynamics (for the same time series of colony growth the first induction occurs irrespective of the value of η), while the shapes of the induction histograms are not significantly modified.

Statistical analysis.
To understand how the induction process is progressing on a statistical level we first address the conditions for the first-induction event. We quantify this by the mean first-induction time-i.e., the time that the first cell switches to the induced state-analogous to the concept of the mean first-passage time for the crossing of a prescribed value of any stochastic process. We therefore use the more common acronym MFPT. As for any given realisation the first-induction time is independent of the production increase η of AI molecules, the relevant parameters are the degree of spatial inhomogeneity and the inverse signal-range parameter. We will not study the effect of variations of the number n c of seed cells, these are assumed to be given by the fixed parameters ρ = 0.0025/μm 2 and R = 100 μm of the Poisson distribution with mean and variance (5). www.nature.com/scientificreports www.nature.com/scientificreports/ We start with the variation of the inverse signal-range parameter α as shown in Fig. 9. For three different values of the standard deviation σ str of the Gaussian daughter-cell displacement distribution, in panel (a) the mean first-induction time is seen to rise with decreasing signal range. This is quite intuitive: the further the signal penetrates the colony away from the producing cell, the more it elevates the cumulative AI concentration due to the lower degree of degradation (for the radial decay of the single-cell AI concentration profile as function of α, see Fig. 2). For a different spread of daughter cells, we see that for higher heterogeneity of the growing colony (smaller value of σ str ) the MFPT is markedly smaller, consistent with our previous argument that in a more clustered colony denser cell regions occur, yielding higher AI concentrations.
For the mean AI concentration 〈 〉 c the decay with α also follows intuition: growing α corresponds to higher AI degradation. While for different colony heterogeneity 〈 〉 c essentially converges in the limit of small α (extended spreading and thus higher insensitivity to disorder) we see that the smallest degree of disorder (σ str = 100 μm) has a slightly (almost within the error bars) higher AI concentration. This is due to the significantly longer MFPT to first induction, see Fig. 9(a), such that more cells have been produced, over-compensating the dilution effect of the larger colony spread. Note that the mean concentration is a relatively bad measure to distinguish the induction threshold for different degrees of heterogeneity σ str , as shown in Fig. 9(b).
The presented data corresponds to the product of many different types of colony clusters and thus includes not only local but also cluster-cluster effects of QS, as investigated in 47 : that is, every colony here can be seen as a superposition of growing clusters of cell lineages that interact with each other in the induction dynamics. Therefore, it is useful to investigate a model colony with spatially isolated seed cells to gain further insight into the degree of locality of the induction process. For this purpose 100 cells were distributed on a grid with a distance of 200 μm to each lattice neighbour. This distance is large enough so that even for very high signalling ranges (very small α) effectively no AI molecules arrive at another cell family. Hence, different cell families will develop independently from each other, until quite mature stages of the colony.
In addition to the MFPT and 〈 〉 c of the first-induced cluster this simulation also delivers data on the mean span 〈r col 〉, the largest distance of any cell in the cluster to the centre of mass of the cluster at the time of first induction as well as its mean number of cells 〈n col 〉 (cluster size). The general trends of the colony with separated cell clusters shown in Fig. 10 are similar to those of the original system displayed in Fig. 9. On closer inspection we see that both the mean first-induction time and the mean AI concentration at the point of first induction have larger values. This fact already reflects the partial interference between clusters. The general behaviour of MFPT and 〈 〉 c in Fig. 10 is not surprising: for short signal ranges (large α or high degradation) the typical concentration is lower and induction on average occurs after longer times. With growing displacement σ str the mean induction time increases, as typical distances between cells are larger. Concurrently, the mean span of the cluster grows with increasing σ str as well as the mean cluster size. The AI concentration of the cluster containing the first induced cell is decreasing with larger σ str as this quantity is only measured per cluster and this scheme disregards the empty space between clusters.
While examining family-separated clusters gives a good impression of the locality of the induction process, the properties of the mean first-induction point do not provide any insight into the later induction process, when additional AI production from already induced cells becomes relevant and the spatial arrangement of the cells starts playing an even more decisive role. To further understand what shapes the induction of new cells it is useful to define a quantity measuring the fraction of induced cells for a statistical population in addition to the fraction χ(t) of induced cells (Fig. S.3). For this purpose we consider the "order parameter" 〈IT χ 〉 defined as the expected time when a fixed χ value (in the interval (0, 1)) is reached. For instance, 〈IT 0.5 〉 is the expected time when 50% of the entire cell population has been induced, potentially an interesting quantity to be studied in experiments. In comparison to the MFPT, 〈IT χ 〉 is a function of the additional production rate η of AI by already-induced cells. We present results for the case χ = 50% in Fig. 11, further plots are gathered in the Supplement. www.nature.com/scientificreports www.nature.com/scientificreports/ We first focus on the case of no surplus AI production (η = 0). As function of the inverse signal-range parameter the general behaviour of 〈IT 0.5 〉 in Fig. S.5 is similar to that of the MFPT in Fig. 9 but with slightly lower standard deviations. Varying the spatial homogeneity parameter σ str similarly shows that more homogeneous colonies (larger σ str values) take longer to reach a certain percentage of induced cells. A clearer distinction can be seen when observing the corresponding mean AI concentration in Fig. S.5(b). While in Fig. 9(b) the behaviour of 〈 〉 c shows a distinct downward trend in dependence on α, when higher fractions of the cells are induced the mean concentration becomes increasingly independent of α. This gradual disappearance of the α dependence can be observed when looking into the same kind of data for lower values of χ (Fig. S.6), as seen for χ = 0.25 in Fig. S.7(a). Here, the effect is weaker than for the case of the MFPT but still more significant than for the case 〈IT 0.5 〉. When larger fractions of a colony are induced and the distances between already induced cells are considerably smaller, degradation effects become less pronounced. In more advanced colonies with higher values of χ another noticeable effect is taking place. While for the point of first induction (see Fig. 9(b)) the mean AI concentration hardly depends on the degree of spatial inhomogeneity set by σ str , we now see a clear separation of the 〈 〉 c values for different σ str , especially for the case χ = 0.75 shown in Figs S.8 and S.9. Higher σ str leads to larger spatial extensions of the colony and thus the respective mean AI concentration is lower. Naively, more homogenous colonies would be expected to have lower mean concentrations because of their larger extension. While first induction is mostly governed by local anomalies of the cell distribution, the spatial-homogeneity parameter gains more influence on the mean concentration at later stages of colony development. Note also that we collected the data with the step of χ being 0.01, but present here the results only for χ = 0.25, 0.5, and 0.75.
Let us now address the influence of the AI production rate increase η. In Fig. S.5(c) as expected, a lower η leads to longer induction times 〈IT χ 〉 while the general behaviour as function of α remains unaffected. The shift of the curves from η = 0 to η = 1 is larger than that from η = 1 to η = 2, pointing at a nonlinear relation between induction time and production increase η, as confirmed for the larger variation shown in Fig. 11. The associated mean AI concentration in Fig. S.5(d) shows clearly lower values for lower levels of additional AI production. For a lower value of χ in Fig. S.7(b) this trend is also observable. The standard deviations of 〈 〉 c are rising with η because The colour coding for different σ str is the same as in Fig. 9. Each data point is taken from 20 simulations runs performed at η = 0.
higher additional AI production overstates the effects of denser regions on the mean AI concentration in a similar way as a longer time evolution would. The detailed analysis of the induction time and mean concentration as function of η in Fig. 11 has a clear message: while the mean concentration increases almost linearly with η, as expected, the decrease of the induction time flattens for larger η values. Thus, increasing the AI amount in the positive feedback loop is mainly beneficial for induction at smaller η values. Too large η values, corresponding to costly cellular processes such as protein production, do not pay off. The influence of the spatial disorder, as shown in Fig. 11(a) in fact is more important than that of the production increase, similar to the effect of α. In that sense physical boundary conditions cannot be compensated by the involved biochemistry.
The associated mean AI concentrations 〈 〉 c shown in Fig. 11(b,d) increase approximately linearly with η. Spatially heterogeneous colonies realised at smaller values of σ str show slightly larger values and standard deviations in Fig. 11(b). This somewhat greater variety of 〈 〉 c can be explained by the number of possible invariant spatial arrangements a system can obtain in regimes of a certain spatial homogeneity. This number is arguably larger in highly heterogeneous systems. Since denser populated regions are more likely to spawn even more new cells, every specific spatial arrangement of cells gets overstated at longer times or larger η. This argumentation is supported when analysing the data for lower values of χ in Fig. S.7(c,d) where the standard deviations are proportionally smaller; and in later time evolution in Fig. S.9(c,d) where the standard deviations are proportionally larger.
In the earlier induction process for χ = 0.25 in Fig. S.7(c,d) the rise of the mean AI concentration 〈 〉 c does not start right above η = 0. Instead, for all levels of spatial homogeneity and all signal ranges presented here it takes a greater amount of additional AI production to influence 〈 〉 c . During the later time evolution of the colony with larger χ values, the rise of 〈 〉 c sets on already at lower values of η, see Fig. S.9(c,d). Colonies with higher levels of heterogeneity show slightly larger concentrations in the earlier process (though there is no clear separation given the resulting standard deviations).
The data in Fig. 11(c) where σ str is held constant at an intermediate value and α is varied from curve to curve demonstrate that the induction time behaves similarly to Fig. 11(a), and it is again shifted by the variation of α but www.nature.com/scientificreports www.nature.com/scientificreports/ not influenced in its general functional form. The associated 〈 〉 c grow with η as well, see Fig. 11(b). We note that for χ = 0.25, as seen in Fig. S.7(d), colonies with large signal ranges α 1/ display slightly higher levels of the mean AI concentration than colonies with small signal ranges. However, this phenomenon vanishes in higher developed colonies, see Fig. S.9(d). The effect inverts during the colony evolution: Fig. S.9(d) displays lower levels of AI concentration at small inverse signal-range parameters α, albeit the standard deviations for α = 0.01/μm 2 are the largest for all χ (due to a higher sensitivitsy of 〈 〉 c on the spatial arrangement of the cells in high signal-range regimes). For smaller α colonies have significantly shorter induction times, that is, less cells have accumulated in the colony when a certain fraction got induced.
The same argument can be used for the variation of σ str . High signal range colonies are more likely to develop a greater number of independent hot spots earlier in the induction process and therefore reach higher mean concentrations earlier. This advantage, however, vanishes over time when lower signal range colonies catch up due to sheer numbers of cells and a higher degree of overlap. This is one of the central results of the current analysis.
The final remaining analysis is the detailed variation of σ str with varying values of α and η. Holding the signal-range constant and using different values of η leads to the results displayed in Fig. S.10. Naturally higher η lead to earlier induction of many cells. The nonlinear effect of increasing η is confirmed in this data series. The time 〈IT χ 〉 grows with the spatial homogeneity parameter in Fig. S.10(a). This behaviour can be seen as well for χ = 0.25 in Fig Phase transition. Since the production of AIs is positively up-regulated after a threshold amount is reached, the vicinity of an induced cell is more likely to get induced, as well, due to elevated concentration of signalling molecules. As seen in the above discussion burst induction is typically a local phenomenon. To emphasise this point we finally study the process of burst induction for a system with a single seed cell. The central order parameter to classify the burst behaviour is the mean relative size 〈ζ〉 of a burst divided by the entire number of cells in the emerging colony at the time of first induction. The main effect we study here is the influence of the up-regulation quantified by η. The resulting data are shown in Fig. 12. We also calculated the expected number of cells in the burst, see Fig. S.11. Figure 12 shows that while for η = 0 hardly more than a few cells are comprised in a burst, with increasing η it becomes more and more likely that an induction burst will comprise significantly larger fractions of the cluster, until the single burst becomes almost global, inducing a large fraction of cells in the cluster. This quite sharp behaviour is reminiscent of a phase transition (or percolation): it is most pronounced in Fig. 12(a,b) for smaller spreads of the daughter cell displacements as well as for the smallest inverse signal-range parameter (longest signal range), in agreement with intuition. The data of Fig. 12 also show that spatially heterogeneous colonies start their phase transition at lower values of η and that their upper limit 〈ζ〉 is closer to 1 (or even exactly one when the signal range is sufficiently high). A larger signal range-corresponding to smaller values of α-always leads to a higher upper limit for all levels of spatial homogeneity.

Discussion
QS enables simple organisms such as bacteria to sample information on the density of same-species or different-species organisms in their local environment. QS controls collective actions of bacterial colonies via production, sensing of, and response to AI molecules. This establishes cell-cell communication and orchestrates collective/communal behaviors of bacterial colonies (sociobiology) 10,13 . QS cascades and reactions at heterogeneous conditions of bacterial biofilms 18 emerge in terms of bacterial density, distribution of nutrients 9 , as well as AI diffusion characteristics. The dynamics of QS under constraints in natural biofilms as well as with fluid flows is important to consider on a model level. Constant flows can wash AIs away thus repressing QS. These and other "external" factors enable genetically identical bacteria to respond to signals and operate individually in different space locations in the colony and at different times. Certain cavities, and crevices of the bacteria-growth medium also impact QS 18 .
We studied the influence of the degree of spatial disorder on the induction dynamics of a growing bacteria colony, as a prototype model for early biofilm formation based on QS when the bacteria mainly form a monolayer on a surface. To this end we set up a Monte Carlo simulations scheme and studied several core quantities such as the shapes of the emerging colonies, their bursty induction behaviour, and the mean AI concentration at the point of induction. We investigated these effects as function of several core parameters: the inverse signal-range parameter (2019) 9:12077 | https://doi.org/10.1038/s41598-019-48525-2 www.nature.com/scientificreports www.nature.com/scientificreports/ measuring the decay of a single cell AI concentration field, the production increase η of AIs upon induction, and the distance by which a daughter cell is displaced from the site of the mother cell. We find that spatial disorder is a main factor in colony growth and induction dynamics. Our results, therefore, add to the current discussion of heterogeneity effects in QS colonies 77,78 .
Our main findings are the highly local induction of clusters in scenarios when an initial colony with random cell placement grows under conditions of short signal range-that is, when the AI concentration field emitted by a given cell falls off over the range of few μm-and relatively local placement of the daughter cells. In this scenario the concentration field becomes highly heterogeneous and the local AI concentration in a given cluster grows much quicker than in the corresponding homogeneous scenario: consequently, the first induction occurs considerably earlier, a new effect observed in this study. We showed that for the studied model parameters this scenario is robust even when the seed cells are evenly spaced: the mixing due to daughter cell placement leads to a disordered geometry of the colony.
A characteristic feature of the induction dynamics is the occurrence of bursts, the simultaneous induction of a number of cells in a neighbourhood. The statistic of bursts turns out to depend crucially on the degree of spatial disorder and the signal range given by the AI degradation. Studied as function of the additional AI production rate, the burst size may show a clear, relatively sharp phase transition from very small bursts to bursts with sizes of the order of the colony size.
Our analysis introduced a family tree such that each cell carries the label of the specific seed cell it stems from. This way simulations can be used to study the degree of overlap between different lineages and effects of the exact growth and induction on a family level. Similarly, we introduced a colour labelling of clusters induced by successive induction cascades. In experiments, it should be possible to endow different cell lines with different colours or other labels. While this may not be possible for large seed colonies in the hundreds, it should be possible to distinguish, say, of the order of ten families in the course of the biofilm formation.
The progress of the colony induction was monitored using different quantities. In addition to the AI concentration and the first-induced cell cluster we introduced the fraction of induced cells at a given snapshot of the colony, as well as the mean first-induction time and the typical times when a given fraction of the colony is induced. Moreover, we studied the distribution of induction levels as function of time, showing distinct features depending www.nature.com/scientificreports www.nature.com/scientificreports/ on the physical parameters of the system. In particular, for given parameters it may happen that while highly heterogeneous colonies are induced much earlier than homogeneous ones, to induce the entire colony it may take longer for the heterogeneous case as major voids in cell distribution have to be overcome upon AI-signal propagation.
An interesting question for biological systems is the economy of resources. Thus, gene expression and molecular signalling processes in cells are typically performed by passive diffusion 79 , as active transport would easily exceed the resources of a cell. The often extremely low concentrations of signalling molecules may be compensated by spatial proximity (geometry-control) 80 between the loci where the molecules are produced and where they are supposed to bind [81][82][83][84] . For the QS scenario we here showed that excessive production of AIs may not lead to more reliable signalling and, thus, allows a similar housekeeping economy for the partaking cells. While the positive feedback loop leading to an up-regulation in the AI production in an induced cell clearly increased the induction dynamics of a growing colony, for even higher up-regulation levels the gain decreases significantly. How much this statement varies with different scenarios of spatial disorder remains to be seen.
The effects of spatial disorder studied here are relatively moderate, taking typical daughter cell displacements of few tens of μm into account, but may provide decisive advantages for cell types with more optimised QS strategies. For mobile daughter cells two scenarios might be considered: (i) cells may actively seek for either already dense or empty areas, thus emphasising or moderating the effects studies here. (ii) Or, cells may choose different, random placement such as long-tailed distributions as used in the location of sparse resources by many animals as well as bacteria 85,86 . Note that even when cells move in a highly randomised fashion, effects such as bulk mediation may lead to non-Gaussian surface displacements at intermediate scales 87 . It should be interesting to see how such scenarios balance the induction versus the spreading behaviour of colony formation by QS. Ultimately, such questions are central for the study of locality versus globality effects in molecular signalling between bacterial cells 88 and in similar systems.
Additional important aspects to study include different, potentially anisotropic dispersal, proliferation, and motility strategies of bacterial cells 8 , non-Fickian or anomalous diffusion of AI molecules 84 , as well as the diffusivity of AI molecules depending on local bacteria density (transport channels and porosity of a biofilm) 8,61,63 .
Other biologically-relevant aspects should include stimulations of cell mortality (when local cell densities become too high to allow flows of nutrients and oxygen), disposal of poisonous metabolic/waste products, of the history of experiencing cell environments 36 , mimicking heterogeneous micro-scale structure of bacterial biofilms/ colonies 9,10 , and QS mechanisms for multiple species of QS-signalling molecules ("interference") 14 . Modelling environment-induced changes in the gene-expression levels and formation of bacterial subpopulations ("phenotypical heterogeneity") 9 are also important to consider.
From the modelling perspective, the effects of excluded volume in dense bacterial colonies will affect diffusion of small molecules through voids and spaces between cells. Naturally, the growth preferences of a bacterial colony from a single seed (provided by a daughter cell) are expected to depend on cell density around the point of first-cell placement as well as on availability of nutrients there. The behaviour of a more complicated model under the conditions of limited amount of nutrients, crowding-dependent growth dynamics and division rates of bacterial cells, and with crowding-dependent diffusion of QS-inducing AI molecules will be interesting to examine. For instance, a constant diffusivity of AI molecules was used in the current analysis, but in denser assemblies of bacteria the cell density will have an effect on D. Our studies of tracer diffusion in crowded media with non-inert obstacles 89 and non-homogeneous crowding 90 have already demonstrated that the effects on free-space tracer diffusivity can be rather dramatic.
Moreover, competition for space, oxygen (aerobic and anaerobic metabolism), and nutrients 10,14,49,70,91 between different cell species with different traits (such as size, spreading properties, or cell-division rates) will give rise to a rich behaviour, for instance, "rock-paper-scissor" competition scenarios [39][40][41] , including inter-species "cross-talk" 50 . We refer to [9][10][11]14 for the discussion of certain evolutionary aspects of QS (relative fitness benefits), selection strategies, and competition-versus-cooperation strategies for public goods in multispecies bacterial colonies. We also refer to 92 for the mechanisms of cooperation, adaptation, and conflict in microbial biofilms.
Naturally, a certain degree of stochasticity of the model parameters is going to quantitatively affect the outcomes and predictions of the current modelling. Importantly, for the cell-division strategies different from the one with a constant rate of division ("timer")-such as the "adder" or "sizer" scenarios 93-95 -the spreading dynamics of a bacterial colony as well as the features of the induction process and general collective behaviour of cells are expected to be affected as well.
While the modelling approach based on individual cells pursued here is very flexible and can accommodate many different features, approaches based on differential equations for chemical concentrations have the advantage that, at least in some limits, analytical solutions can be obtained 26,47 . In addition to these two strategies network-based approaches have been suggested 96 . We finally note that some of the results observed here may have relevance in morphogen-gradient based problems, in which the spatial distribution of the concentration of stimulating molecules can trigger a certain response of tissues.

Data Availability
The data sets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.