Spatial organization has a strong influence on the development and dynamics of biological systems (Kreft et al., 1998; Lardon et al., 2011; Halsted et al., 2016; Hansen et al., 2016; McNally et al., 2017; Timm et al., 2017). The factors leading to changes in organization of multicellular assemblages at fine spatial scales are not well characterized or understood, however, they are expected to guide the succession of community development and function toward specific outcomes (Liu et al., 2009; Cline and Zak, 2015; Dini-Andreote et al., 2015). The organization of distinct microbial populations can be shaped by physical and chemical processes, and affect important activities such as antibiotic resistance, efficient energy conversion, C and N cycling and quorum sensing (Ginovart et al., 2005; Gras et al., 2010, 2011b; Sahari et al., 2014; Wang and Ma, 2014; Koonin and Wolf, 2015; Biteen et al., 2016). Microbe-microbe interactions can also depend on direct and indirect competition for resources between different community members (Kreft, 2004; Hellweger et al., 2008; Borenstein et al., 2015; McNally et al., 2017). The microscale/local transport of essential microbe-derived metabolites and cell-to-cell competition are likely to be strongly influenced by spatial confinement and individual cell behavior in the environment (Lardon et al., 2011; Pintelon et al., 2012; Vogel et al., 2015; McNally et al., 2017). Consequently, investigating the complexity of these processes and emergence of unique behaviors requires the combination of experimental and computational tools that can be used to explore the impact of spatial organization, while correlating individual microbial behavior and interactions to specific outcomes (Dini-Andreote et al., 2015; Zhu et al., 2015; Hansen et al., 2016).

Cells can compete directly with surrounding species through physical contact, and in more specialized cases, are capable of transferring toxic effector proteins to susceptible cells. The Type VI secretion system (T6SS) is an important example of such a pathway, being responsible for the assembly of a pilus apparatus that can be used to contact neighboring cells and potentially induce cell death (Hood et al., 2010; Chou et al., 2012; LeRoux et al., 2012). Hood et al. (2010) showed that the H1-T6SS of Pseudomonas aeruginosa is required to direct the injection of toxins from T6SS active cells (T6SS+) into T6SS-susceptible cells (T6SS−) that lack immunity. Other important secretion systems such as H2- and H3-T6SS in P. aeruginosa direct toxins preferentially to eukaryotic cells. However, because the H1-T6SS toxin is preferentially directed toward other bacteria, it is particularly well suited for studies of contact-mediated interactions between neighboring and competing prokaryotes (Mougous et al., 2006; Sana et al., 2016, 2017). T6SS interactions in mixed microbial populations also play an important role in the regulation of more complex biological processes and microbial community dynamics (Russell et al., 2014; Verster et al., 2017). For instance, the T6SS interactions occurring amongst commensal bacteria in the mammalian gut microbiome have been shown to modulate community composition and interactions, as well as provide a mechanism for defending commensal bacteria from invading pathogens (Hecht et al., 2016). Furthermore, these T6SS interactions are highly active and prevalent, where > 10⁹ T6SS firing events (i.e., predicted pilus injections) min⁻¹ g⁻¹ colonic contents can occur and nearly 25% of human gut microbiota have been shown to encode a T6SS pathway (Wexler et al., 2016; Sana et al., 2017).

Using two-member communities as a model system of T6SS interactions in the laboratory, Borenstein et al. (2015) demonstrated that established colonies of T6SS− Escherichia coli could survive contact with T6SS+ Vibrio cholerae. Agent-based modeling (ABM) simulations further showed that T6SS− cells could survive T6SS+ attack when placed in situations of nutrient limitation and relatively slow growth rates, and could even outcompete the T6SS+ cells, as long as T6SS− cells were able to establish microcolonies within the mixed community (Borenstein et al., 2015). These results demonstrate the importance of spatial confinement and local organization on cell growth and survival. Thus, competition between neighboring microbial cells and spatial confinement are expected to drive changes in cell assemblage and organization (Borenstein et al., 2015; Halsted et al., 2016; Hansen et al., 2016).

Numerous advances in our understanding of cell-to-cell behavior and interactions at fine spatial scales have stemmed from the use and development of nano/micro-fabricated platforms (Wang et al., 2013; Yamazaki et al., 2014; Swennenhuis et al., 2015; Xue et al., 2015; Hansen et al., 2016; Zhang et al., 2016; Timm et al., 2017; Yeh et al., 2017). Timm et al. (2017) used a microwell array platform to study the contact-mediated T6SS interactions of P. aeruginosa. The microwell array platform enabled high-resolution and high-throughput imaging of mixed T6SS+ and T6SS− cells growing under spatial confinement within microwells, with well diameters ranging from 20 to 100 μm and 5 μm depth. Interpreting the results of these cell-to-cell interactions with simplified analytical models of overall growth within each well becomes challenging and potentially unreliable when trying to capture the complex interactions reflected by spatial organization of microorganisms within the microwells. Alternatively, ABM simulations can capture how changes at the level of individual microbial interactions lead to changes observed at the community and microcolony levels. In conjunction with laboratory experiments ABM simulations can be used to infer and test important growth parameters that impact spatial organization within colonies (Borenstein et al., 2015).

In this study, we have developed an ABM model around experimental data obtained from a microwell array platform. We use the model to interpret spatial organization patterns of P. aeruginosa mutants growing under spatial confinement. The novelty of our approach relies on the high throughput nature of both the experiment and ABM simulations, which allows investigating how the initial ratio of community member abundances, initial growth location and T6SS interactions affect spatial organization during growth. The model is written in the language NetLogo (Wilensky(1999), NetLogo, http://ccl.northwestern.edu/netlogo/; Center for Connected Learning and Computer-Based Modeling, Northwestern University, Evanston, IL) and is linked to a computational framework that permits submitting many calculations in parallel for different initial parameters, where each combination of parameters can be conceptualized to represent a micro-environment of interest. The ABM model has been deployed in the Compute and Data Environment for Science, CADES (http://cades.ornl.gov/), which also stores the relevant experimental data used during fitting routines. We find that key parameters, such as spatial constraints, local concentrations of growth substrates/metabolites and associated rate constants alter the impact of P. aeruginosa Type VI secretion activity on the spatial organization of cells in confined environments.

Previously, Timm et al. (2017) demonstrated that the GFP ΔretS mutant T6SS+ and the m-Cherry ΔretS/Δtse/i1-6 deletion mutant T6SS− that is susceptible to Type VI secretion, formed discrete microcolony assemblages during co-growth in microwells. Observed well populations were heterogeneous (Figure 1A) with distinct assemblages of individual species forming across the wells. Instances of T6SS− microcolony formation were unexpected based on the hypothesis that T6SS+ cells would dominate each well environment due to their directed injection of toxic effector proteins into susceptible cells. In the present study, we used ABM simulations (Figure 1B) to examine how spatial confinement and Type VI secretion can lead to growth conditions that allow the formation of microcolony assemblages and enable susceptible T6SS− cells to persist.

In Figure 3, we show how the experimental values for a′, μ′, τ change with the initial T6SS+: T6SS− seeding ratio for both types of cells in 30 μm diameter wells. Although the initial ratio of cells in the solution used to seed the wells is 1:2, natural variability in the number of cells of each type that seed into the wells provides initial ratios ranging between nearly 0 (all T6SS−) to 3:1 (3 times more T6SS+ than T6SS−). For both types of cells, τ is practically independent of the initial ratio, whereas a′ and μ′ for T6SS+ linearly decrease as the T6SS+:T6SS− ratio increases. Likewise, μ′ and a′ values decrease as cultures have an increasing number of T6SS− cells compared to T6SS+. For both T6SS− and T6SS+ the maximum number of cells per initial number of cells (a′), as well as the maximum rate of increase per initial numbers of cells (μ′), decreases as the corresponding cell type is initially present in excess. As seen in Figure 4, for each 30 μm well, the initial density is relatively constant, independent of the initial T6SS+:T6SS− ratio. Cells are, for the most part, segregated into domains that contain either T6SS + or − cells. The latter is observed both in experiment and simulations (Figures 1A,B, respectively), though some overlap can be seen within the experimental images (Figure 1A, yellow region). We speculate that larger clusters or colonies of cells will grow more slowly on a per cell basis because of nutrient limitations that may occur at the center of those microcolonies relative to growth at the edges. Thus, even though the overall density within the wells is the same, local microcolony size may influence the maximum growth rate per initial number of cells (μ′) and maximum number of cells per initial number of cells (a′). A more detailed image analysis strategy that allows quantitative description of microcolony size and patchiness within wells throughout these experiments is warranted and under development.

During growth on the experimental microwell platform, both T6SS + and − cells showed a general preference for seeding at the well edges regardless of well size (Timm et al., 2017). During ABM-simulated growth in a circular well with preferential seeding at the boundaries, microcolony formation qualitatively mirrors the experimental results of Timm et al. (2017) (Figure 1) though, again, a more quantitative analysis of microcolony size and patchiness within wells is needed to make a direct comparison. In the simulations, the cell placement is biased so that there is a higher probability of cells beginning the simulation at the edges of a well. Cells seeded near the edge of wells, in both the experiments and the simulations, propagate toward the interior because of confinement imposed by the well edge. Based on the contact mediated pathogenesis associated with T6SS interactions, we expect that subsequent cell-to-cell interactions would help maintain bacterial domain segregation and minimize co-localization.

In Figure 5 we have plotted the slope of the straight lines, i.e., a′/ratio, μ′/ratio and τ/ratio, for each well size vs. the well diameter. In the case of T6SS+ cells, the slope of these lines shows more variation than that seen for T6SS− cells. The reason for this remains unclear and may result from T6SS+ cells growing outside the wells and outside the analytical region of interest (see Supplementary Movie 1). Also, fluctuations in GFP expression or loss of GFP intensity may contribute to this variation across experiments. For T6SS− cells, the sensitivity (slope, μ′/ratio and a′/ratio) of maximum growth rate per initial number of T6SS-cells and maximum number of cells per initial number of cells, increases slightly with well diameter. Figure 4 shows that the initial seeded cell density per well on the experimental platform decreased as a function of increasing well diameter between 15 and 100 μm. One might expect that at lower cell densities, the sensitivity to initial ratio of T6SS+ to T6SS− cells would decrease because lower densities should correspond to fewer T6SS+ to T6SS− interactions. However, if increasing microcolony size is suspected to cause reductions in μ′ and a′, as described above, larger wells could facilitate formation of these larger domains by (i) allowing unimpeded growth of microcolonies and (ii) increasing the possibility for the seeding of larger microbial aggregates or flocs from solution. Again, more detailed image analysis techniques should facilitate future investigations using the array platform.

Using the ABM model, we performed a similar study of microcolony formation in a population of T6SS+ and T6SS− cells, where the aggressiveness of the former was changed from high to moderate to low. For simulating the experimental conditions closely, we computed the experimental distribution of T6SS+:T6SS− initial ratios and densities and used similar initial values in the simulations. The results of the simulations are shown in Figure 6 for a well size of 30 μm, where we have plotted a′, μ′ and τ for initial T6SS+:T6SS− ratios that resemble those obtained experimentally. When the aggressiveness was 1 (i.e., each contact occurring between a T6SS + and a T6SS - cell results in the death of a T6SS− cell), the behavior of a′, μ′ and τ for T6SS+ resembles that seen experimentally (compare Figure 6A to Figure 3A). The opposite is true for T6SS− cells (Figure 6B vs. Figure 3B), and the results indicate that the fewer T6SS− there are present, the more poorly T6SS− cells grow. The trend seen in Figure 6B for T6SS− cells can be explained as follows: while T6SS+ cells are not affected by their own aggressiveness, T6SS− cells are, and although it is theoretically advantageous for T6SS− cells to have relatively fewer neighboring T6SS− competitors, this advantage is offset by a high rate of lysis caused by the highly aggressive T6SS+ cells. Lowering the aggressiveness of T6SS+ from 100 to 10% does not significantly change the behavior of a′, μ′ and τ for both T6SS+ and T6SS− mutants (Figures 6C,D). However, if the aggressiveness is lowered to 1% (Figures 6E,F), then a′ and μ′ of T6SS− cells follow the same trend as in the experiment (Figure 3B). At this low killing rate, however, the growth rate of T6SS+ cells decreases (compare Figure 6E to Figure 3A). Specifically, as the aggressiveness of T6SS+ cells is reduced, the linear fit of a′ intersects the y-axis at a lower point. Thus, at lower abundances and lower levels of aggressiveness, T6SS+ cells can no longer effectively compete with T6SS− cells.

The T6SS of P. aeruginosa is an important biological model for understanding how cell-to-cell contact directs the succession and organization of microbial communities (Robinson et al., 2009; Hood et al., 2010; Sarris and Scoulica, 2011; LeRoux et al., 2012; Das et al., 2013). As mentioned in the introduction, T6SS interactions play a significant role in the regulation of microbiomes, which has important implications for biomedical and pathogen research, particularly for understanding mammalian gut microbiomes, and also environmental biogeochemistry relevant to native microbial interactions with plants and soil. However, the factors leading to changes in organization of microbial cells at fine spatial scales, driven by T6SS interactions, are not well characterized or understood. Recent laboratory investigations and ABM simulations indicated that established T6SS− colonies of Escherichia coli could persist during cell-to-cell interactions with Vibrio cholerae T6SS+ cells (Borenstein et al., 2015). The results of Timm et al. (2017) further suggested that spatial confinement, as well as T6SS activity between growing effector and susceptible P. aeruginosa mutants, could potentially direct cell organization in micro-colonies and affect the survival of susceptible cells. In the present study, building upon additional analysis of the complete dataset of Timm et al. (2017), and in combination with ABM simulations, we provide supporting evidence that both spatial confinement and T6SS activity can lead to changes in the organization and persistence of P. aeruginosa.

We found that discrete zones of clearing occurred around T6SS− cell assemblages during co-growth with T6SS+ cells in ABM simulations (Figure 1; see also Supplementary Movie 2). This cell-to-cell organization of T6SS− cells, surrounded by a zone of clearing, is consistent with T6SS-induced cell lysis at the boundary between both P. aeruginosa mutants (Hood et al., 2010; Borenstein et al., 2015). This zone of clearing provides a mechanism of P. aeruginosa cellular organization, as previously observed in Si-based microwell arrays (Timm et al., 2017). We speculate that during growth of both mutant strains, these buffer zones can occur randomly during growth, perhaps forming safe-pockets for susceptible cells to continue growing, and can become more defined as microcolonies of both species expand and interact at their outer boundaries. Fitting of the complete experimental dataset indicated that starting at the apparent peak in cell growth for both strains, a general decay in T6SS+ GFP signal intensities began, while T6SS− m-Cherry intensities subsequently remained more persistent over time (Figure 2). Borenstein et al. (2015) demonstrated that more-established microcolonies of T6SS-susceptible cells can potentially survive T6SS attack, which helps explain the persistence of susceptible P. aeruginosa mutants as deduced from the fluorescence intensities taken from our experimental data. We also observed that T6SS+ cells could outgrow a well once the interior of the well cavity had become nearly filled by growing cells (Supplementary Movie 1); this may explain, to some extent, the sharp decay phase generally observed for T6SS+ GFP intensities (Figure 2).

We found that the initial seeded cell density per well on the experimental platform decreased as a function of increasing well diameter between 15 and 100 μm (Figure 4), but cell density did not have an apparent effect on cell organization during different growth simulations, which is consistent with the results of Timm et al. (2017) that demonstrated microcolony formation across all well sizes between 20 and 100 μm diameters. The correlation between initial cell density after seeding and well size likely reflects the preparation of the experimental microwell platform. For instance, following the experimental cell seeding step (Timm et al., 2017): (1) slight drying of the aqueous culture media before contact with the nutrient agarose cover; (2) difficulty rinsing cells from smaller diameter wells during the final water rinse step; or (3) a larger side-wall to floor area ratio per well could have affected initial cell densities such that smaller wells were more densely packed than larger wells, particularly at well edge boundaries. However, qualitatively, we found that spatial organization into distinct T6SS mutant assemblages during experimental and simulated growth was not strongly influenced by cell density or close packing. Future quantitative analysis of assemblage size and spacing for different well sizes may reveal a more defined mechanism. Densely packed cell assemblages have been shown in previous studies to follow similar biological-phase separation where distinct microcolony formation is favored regardless of cell-to-cell density in spatially confined environments (Tolker-Nielsen and Molin, 2000; Berk et al., 2011; Borenstein et al., 2015; Cutler et al., 2015; McNally et al., 2017).

We generally found the impact of well size to be negligible for the size ranges explored in experiments, see Figures 5A,B. Well size did not have a significant impact on overall growth rates per initial cell number or maximum growth rate per initial cell number. This was unexpected. Indeed, in well diameters < 25 μm, competition for resources and cell-to-cell interactions would have been expected to suppress T6SS− growth. Reductions in spatial confinement within larger wells would, in principle, allow T6SS− cells to grow more efficiently with increasing well diameters, reducing the likelihood of encountering T6SS+ because of the lower seeding densities and more available area. In wells of 45 μm diameter and greater, at much lower initial densities (Figure 4), the individual T6SS− colonies may have had the potential to develop with less competition and become more established before interacting with the more aggressive mutant strain. In this case, the perimeter of T6SS-interactions around a colony would be overshadowed by the more established interior of each mutant colony. In other words, with larger wells above 40 microns, T6SS-killing should have become secondary to the size of mutant colonies by the time they interact at their edges. Clearly, a more systematic study of micro-colony size and distribution is needed to understand these results. The behavior of T6SS+ cells across well sizes displays variation in the data that makes it difficult to draw specific conclusions about T6SS+ growth as a function of well size.

Finally, average a′, μ′ and τ for each well size vs. the entire well size distribution were also calculated, see Figure 7; the results obtained in a mixed population of T6SS+ and T6SS− cells were compared to those obtained in control experiments comprising only one cell type. As seen in Figure 7D, for the control experiments, the average a′, μ′ and τ of T6SS− cells are practically insensitive to the well size, whereas the same data for T6SS+, Figure 7C, shows variation and an increase at larger well sizes. In mixed populations, the data for T6SS− cells shows a similar trend, Figure 7B, although the error bars are larger, illustrating the interactions with T6SS+ cells. The data for T6SS+ cells, Figure 7A, shows even larger error bars, which is surprising because these cells should not be negatively impacted by the presence of T6SS− strains. As mentioned above, inspection of the data reveals that on some occasions T6SS+ cells can outgrow/leave the well boundaries (Supplementary Movie 1). We believe this is one of the primary reasons for the large variation observed in Figures 7A,B. Consequently, whether or not T6SS+ cells outgrow or escape the wells should also affect the growth of T6SS− cells remaining within the same wells.

In this study we have developed an experimental-ABM framework that can be used to interpret unique spatial-organization patterns of P. aeruginosa cells growing under spatial confinement. The ABM model developed here, although qualitative, is capable of showing microcolony formation regardless of initial density, or whether the bacteria prefer to begin growth at the well edges, which is consistent with other recent studies that have examined different microbial species under T6SS interactions. Our model was also capable of reproducing the behavior of a′, μ′ and τ of T6SS+ cells for a particular well size, and the same was true for T6SS− cells once the aggressiveness level of T6SS+ cells was lowered. As such, this model can be used to extract information regarding aggressiveness levels, amount of available resources, and rate of consumption of nutrients, as well as how all these variables affect the growth of the bacterial colony. Yet, there are some uncertainties that the current model does not take into account. Future work will focus on optimizing the model by identifying the most essential growth parameters and developing a more quantitative description of variables used for running the simulations. Further, the ABM model is capable of investigating 5000 wells in 30 min, and in connection with the microfluidic platform, constitutes a powerful framework to connect microbiological experiments to ABM simulations, while improving the ABM models to more accurately reproduce the experimental observations. Finally, this new microfluidic-ABM framework could be used in the future to predict the types of microcolonies that are likely to develop when different microbial species are mixed, which is expected to advance our understanding of microbial ecology at fine spatial scales, as well as mechanistically describe how microbial succession occurs in nature and shapes environments of interest.

All authors were involved in the writing and proofing of the manuscript. JW, AT, MH, and SR were responsible for the design, execution and interpretation of the experimental platform and ABM results. PD, JDA, MG, CP, XP, and MF-C were responsible for implementing the ABM model, doing the simulations and interpreting the results.