In this work we make use of the interaction network of a set of 78 bacterial strains isolated and studied by Pérez-Gutiérrez et al. (2013). The strains in this set were isolated from 5 different samples taken from the superficial sediment of Churince pond in Cuatro Cienegas, Mexico. Since the isolating methodology involved subjecting the samples to thermal shock, all of the isolated strains came out to be thermo-resistant, and most of them belong to the genus Bacillus. All 78 × 78 pairs of bacterial strains were cultured in Petri dishes to test for antagonistic interactions. The resulting antagonism matrix is reported in Pérez-Gutiérrez et al. (2013) and reproduced in Figure 1. In this figure, bacterial strains are organized in decreasing order according to their Aggressiveness Index (number of other strains antagonized by a given strain, minus number of other strains antagonizing it). The ID numbers given in this work to all strains, the labels employed by Pérez-Gutiérrez et al. (2013), and the corresponding aggressiveness indexes are tabulated in Table 1.

In order to simulate the evolution of a bacterial community interacting according to the antagonism matrix reported in Pérez-Gutiérrez et al. (2013), we developed a computational algorithm as follows:

A square grid of 200 × 200 cell represents the pond soil. Each cell in the grid corresponds to a small surface that can be either empty or colonized by at most one bacterial strain.

In order to test how the system dynamics depends on the architecture of the antagonism interaction-network, we constructed alternative antagonism matrices (and employed them to repeat the analysis of the system dynamics) following two different strategies:

Random antagonism matrices were constructed by randomly generating N unidirectional antagonistic interactions between the 78 isolated bacterial strains, with N the number of links in the original matrix. We took care of avoiding self-antagonism.

To investigate whether, as claimed by Pérez-Gutiérrez et al. (2013), antagonistic interactions have a lead role in the spatial distribution of bacterial communities, as well as in the diversity differences found across sites, we took the antagonism matrix reported therein, and used it to implement the algorithm described in the Section 2. Recall that, in such algorithm, a 200 × 200 square grid represents the pond, and each cell in the grid corresponds to a small surface in the pond soil that can either be empty or colonized by at most one strain. Initially, the grid was “inoculated” with 50 individual colonies of every bacterial strain, randomly distributed across the grid (see Figure 2A). That is, the grid was seeded with 78 × 50 = 3900 cells in total. The system was then iteratively evolved in time for 500 round cycles of the algorithm. In these initial simulations, we assumed for all strains (μ) that the probability that the bacterial strain in a grid cell colonizes a neighboring empty cell in an algorithm step is 1 (P_g(μ) = 1), and that the probability that an inhabited grid cell is emptied because the corresponding colony dies is zero (P_d(μ) = 0). To clarify the role of antagonism in the evolution of the system, an aggressiveness index was defined as follows: the aggressiveness index of a given bacterial strain, AI(μ), equals the number of other strains antagonized by it, minus the number of other strains that antagonize it.

In Figure 2 we present the results corresponding to a single simulation. We carried out several simulations and obtained equivalent results in all cases. Namely, an emergent spatial arrangement of bacterial communities in patches spans the grid after a few iteration steps, and eventually a stationary state in which the patches no longer change is reached (see Figure 2B). We note that the formation of patches is consistent with the fact that different strains were isolated from different samples (Pérez-Gutiérrez et al., 2013), despite the habitat homogeneity. Moreover, we can see in this particular simulation, but this is a consistent observation, that patches of vulnerable strains are always surrounded, and shielded from the most aggressive strains, by patches of other strains that are both resistant to the aggressive ones and non-antagonistic toward vulnerable bacteria. This result suggests a mechanism for the survival of low-antagonistic and vulnerable bacteria, but also offers a plausible explanation for the observation that antagonistic interactions within a sampling site are on average less frequent than interactions across sites (Pérez-Gutiérrez et al., 2013).

To have a deeper understanding of the system dynamics we measured patch sizes in the stationary state of 100 different simulations. The results are summarized in Figure 3A. We can appreciate there that a positive correlation exists between the mean patch size and the aggressiveness index of a given strain. That is, more aggressive bacteria tend to form larger patches, whereas the patches of more vulnerable bacteria are consistently smaller. We further measured the total population of every strain at the end of 100 different simulations and computed the corresponding averages. A scatter plot representing mean final population size vs. aggressiveness index is shown in Figure 3B. Observe that total final population is also positively correlated with the aggressiveness index.

Finally, we followed the evolution in time of all the strain populations and the results (averaged over 100 different simulations) are shown in Figure 3C. Observe that, initially, all the inoculated colonies grow steadily, but eventually, a stratification of them according to their AI value becomes apparent. As far as we understand, this happens because reduction of the distance between colony borders increases the frequency of conflicts. After about 200 simulation steps, a stationary regime, in which the population size of all bacteria strains remains constant, is reached.

Since the spatial arrangement of bacteria in stagnant community patches appeared to be necessary for the preservation of bacterial diversity, we carried out an alternative experiment with constant shuffling (mixing) of the grid cells. In these experiments, we carried out simulations as previously described (see Figure 4A), but every 10 simulation steps the grid cells were randomly shuffled (different inter-shuffling times were also considered, but all of them led to equivalent results). We observed that all susceptible bacterial strains (most of them with low antagonistic levels) become extinct, and they are replaced by larger populations of medium- and high-level antagonist strains (see Figures 5B,C). Furthermore, the surviving bacteria do not form patches. Instead, single-grid-cell colonies are randomly distributed (in a well mixed fashion) across the grid (see Figures 4B, 5A). Even though we shuffled the grid at regular times, the first shuffling events have the most notorious effects: they rapidly cause the extinction of the sensitive strains. After this, the populations of surviving bacteria reach a stationary value, independently of the changing spatial distribution (see Figure 5C). A correlation between stationary populations and AI values is still present for the surviving strains (see Figure 5B).

The above results indicate, to our consideration, that periodic shuffling of the grid cells not only precludes patch formation (by randomly relocating individual colonies of all strains) but also ensures that all susceptible strain colonies eventually get in contact with aggressive strains and are thus driven to extinction.

To answer the question of whether the results discussed in previous subsections depend on the architecture of the antagonistic-interaction network (which is determined by the antagonism matrix), we generated random and experimental-like interaction matrices (see Section 2). Experimental-like interaction matrices have the same number of high-level, medium-level, and low-level antagonistic strains, and the distribution of exerted-antagonism and received-antagonism links are very similar to those of the experimentally obtained matrix. On the other hand, random antagonism matrices were built by linking, via antagonism interactions, couples of bacterial strains chosen at random (avoiding self-antagonism). The number of antagonism links in the random and the experimentally-obtained matrices are always the same.

When we repeated the simulations using experimental-like interaction matrices, we were able to recover all the previously described results, with and without shuffling. Nonetheless, the simulations with random interaction matrices rendered quite different results. Community patches can still be observed in the non-shuffling simulations, but the final bacterial population size is less disperse than in the simulations carried out with the experimentally-obtained or the experimental-like antagonism matrices. Furthermore, in the shuffling simulations, most bacterial strains become extinct and no correlation is observed between AI and the stationary population size (see Supplementary Material).

In conclusion, the architecture of the interaction network plays a very important role in the emergence of the community patch structure. In particular, the characteristic that seems to be essential is that bacterial strains can be classified in three different classes: aggressive (highly antagonistic and resistant to other strains), neutral (barely antagonistic and resistant to aggressive strains), and vulnerable (non-antagonistic and sensitive to aggressive strains).

So far, we have assumed in our simulations that all bacterial strains grow at the same rate, whenever they have available space. Recall that in our model the growth rate of strain μ is determined by probability P_g(μ). Moreover, we considered equal initial populations in all of our simulations. In order to have a more complete picture of the studied biological phenomena, we repeated our simulations by taking into consideration distinct growth probabilities and initial populations.

Based on previously reported observations (Bohannan et al., 2002; Kneitel and Chase, 2004; Cadotte, 2007), who have shown that changes in the flow of metabolic energy and in the cellular machinery need to be done by a bacterium in order to become either highly antagonist or fast growing, we assumed that P_g(μ) is inversely proportional to AI(μ) (the strains with the largest and the smallest AI values have P_g = 1 and P_g = 0.5, respectively), and repeated the previously described simulations. According to our results, the new simulations differ from the previous ones in that low-level and medium-level antagonistic strains have now larger stationary populations, at the expense of high-level antagonistic strains (see Supplementary Material).

To account for variable initial populations, we took the final population distributions of the simulations described in Section 3.1, and used them to compute the initial strain populations of a new grid, with the constraint that the most populous strain initially occupies 1% of the grid cells. After running the simulations, the only notorious difference we observed, with respect to those corresponding to constant initial populations, is that having larger initial population has a positive effect on the corresponding stationary population. This happens regardless of whether we consider constant or variable growth rates (see Supplementary Material).

Finally, we considered in our model a death rate related to causes independent of bacterial antagonism. According to Servais et al. (1985); Pace (1988), in aquatic environments there are several events that generate the reduction of certain populations of bacteria, due primarily to: biotic factors (viral infections, starvation, senescence, and allelopathy) and abiotic factors (environmental changes such as radiation, temperature changes and toxic compounds). To take this into account we assumed that, in every step of the algorithm, a grid cell inhabited by strain μ has a probability P_d(μ) of being emptied, independently of growth and antagonism events.

Initially, we considered a constant P_d(μ) = 0.01 for all strains. With this assumption, we were able to recover all the results described in previous sections, with the exception that the most favored strains (those with the largest stationary populations) are those with intermediate AI values. I.e., those strains that neither antagonize many other bacteria nor are susceptible to be antagonized by the most aggressive ones (see Supplementary Material). When we repeated the simulations with P_d(μ) = 0.1, we obtained similar results, except that the strains that had the smallest stationary populations in the previous simulations were driven to extinction (see Supplementary Material). Interestingly, when a non-zero death rate is accounted for, the achieved stationary state is not stagnant any longer. Community patches are formed, but they do not remain the same once the stationary population levels are achieved and the grid is full. On the contrary, they slowly change their size and shape, and move around the grid.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.