# Antagonistic interactions are sufficient to explain self-assemblage of bacterial communities in a homogeneous environment: a computational modeling approach

^{1}Unidad Profesional Interdisciplinaria de Ingenierías Guanajuato, Instituto Politécnico Nacional, Silao, Mexico^{2}Departamento de Ingeniería Genética, Unidad Irapuato, Centro de Investigación y Estudios Avanzados del IPN, Irapuato, Mexico^{3}Unidad Monterrey, Centro de Investigación y Estudios Avanzados del IPN, Apodaca, Mexico

Most of the studies in Ecology have been devoted to analyzing the effects the environment has on individuals, populations, and communities, thus neglecting the effects of biotic interactions on the system dynamics. In the present work we study the structure of bacterial communities in the oligotrophic shallow water system of Churince, Cuatro Cienegas, Mexico. Since the physicochemical conditions of this water system are homogeneous and quite stable in time, it is an excellent candidate to study how biotic factors influence the structure of bacterial communities. In a previous study, the binary antagonistic interactions of 78 bacterial strains, isolated from Churince, were experimentally determined. We employ these data to develop a computer algorithm to simulate growth experiments in a cellular grid representing the pond. Remarkably, in our model, the dynamics of all the simulated bacterial populations is determined solely by antagonistic interactions. Our results indicate that all bacterial strains (even those that are antagonized by many other bacteria) survive in the long term, and that the underlying mechanism is the formation of bacterial community patches. Patches corresponding to less antagonistic and highly susceptible strains are consistently isolated from the highly-antagonistic bacterial colonies by patches of neutral strains. These results concur with the observed features of the bacterial community structure previously reported. Finally, we study how our findings depend on factors like initial population size, differential population growth rates, homogeneous population death rates, and enhanced bacterial diffusion.

## 1. Introduction

Microbial ecosystems have proved to be excellent frameworks to understanding ecological systems (Prosser et al., 2007). New ecological theories have arisen from microbial ecology due to its simplicity, controllability, replicability and experimentally-required times (Jessup et al., 2004). Remarkably, much of that progress has been achieved by means of simplified theoretical models (Momeni et al., 2011). Most of these models account for the interaction of only a few microbial populations. This is an advantage, because simple models can be more easily studied both numerically and analytically, but also a limitation, because such models oversimplify biological reality. Nowadays, thanks to the rapid increase of computer power, it is possible to investigate the dynamics of larger sets of interacting populations (Costello et al., 2012).

Organisms in an ecosystem are affected by abiotic and biotic factors. Abiotic elements refer to non-living chemical and physical factors in the environment, whereas biotic factors are living or once-living organisms in the ecosystem. The influence of abiotic factors on a community spatial-temporal dynamics has been thoroughly studied (Dunson and Travis, 1991). Conversely, biotic factors have been mostly disregarded. Nonetheless, recent studies have evidenced that biotic interaction networks play very important roles (Hooper et al., 2005), specially in terms of biodiversity maintenance and energy flow (Raes and Bork, 2008; Eisenhauer et al., 2012).

Diverse biotic interactions have been experimentally described in nature, and some others have been hypothesized from mathematical modeling results (Evans et al., 2013). At microbial scales, an outstanding panoply of interaction mechanisms has been observed (Prasad et al., 2011), favoring competition over cooperation (Foster and Bell, 2012). Antagonistic interactions (also known as interference competition) have been studied by Czárán et al. (2002), using the simplified model: Killer (K), Resistant (R) and Sensitive (S). This interaction loop implies an associated metabolic cost (and a concomitant growth-rate reduction) for being either an antibiotic producer (K) or a resistant strain (R). Thus, in terms of proliferation, S outcompetes R and R outcompetes K. This model, which is able to sustain bacterial populations, gave rise to new studies aimed at figuring out how microorganisms use and evolve bacteriocin and antagonistic molecules to sustain biodiversity in structured ecosystems (Kerr, 2007).

The shallow water system of Churince, located in Cuatro Cienegas, Mexico, sustains an impressive microorganism biodiversity, mainly due to its geological history (Souza et al., 2006). Even with scarce nutriments in the water, microorganisms have proliferated and evolved to use any available component around them, as has been seen in other oligotrophic environments (Kuznetsov et al., 1979). It is well known that some bacteria, such as *Bacillus spp*., are able to synthesize antagonistic molecules (Abriouel et al., 2011). Cuatro Cienegas' bacteria are not an exception because, as previous studies show, these organisms use a variety of bacteriocin molecules to annihilate competing neighbors (Pérez-Gutiérrez et al., 2013; Aguirre-von Wobeser et al., 2014). All of this, together with the fact that the sediment of this natural setting is static and homogeneous regarding its physicochemical conditions, make it an excellent candidate to study how biotic interactions affect the spatial-temporal arrangement of microbial populations.

In a previous study (Pérez-Gutiérrez et al., 2013), we isolated 78 bacterial strains (most of them from the genus *Bacillus*) from 5 different sampling sites across Churince pond, and tested them for one-to-one antagonistic interactions. Among others, we obtained the following results which are not self evident from the point of view of the most commonly accepted ecological theories:

• Bacterial strains are not homogeneously distributed across the pond, in spite of the pond's physical-chemical conditions being homogeneous and pretty much stationary.

• Antagonistic interactions are more frequent across sampling sites than within them.

In 2013 Pérez-Gutiérrez et al. (2013), we hypothesized that microscopic microbial antagonistic interactions may be responsible for shaping bacterial communities at the macroscopic scale, and that this may suffice to explain the above-enlisted observations. The present work is aimed at proving the feasibility of such hypothesis from a mathematical modeling perspective. To that end we decided to model the pond as a square grid, each of whose cells represents a small area that can be colonized by at most one bacterial strain. The dynamics of the grid cells are then modeled as a set of rules derived from the antagonism matrix experimentally determined in Pérez-Gutiérrez et al. (2013). We decided to employ this modeling strategy because of its adequacy given the available experimental data, and because it has been employed to demonstrate how local interactions can give rise to complicated global patterns, like those observed in natural systems (Gardner, 1970; Hogeweg, 1988; Iwasa et al., 1998; Sarkar, 2000; Wootton, 2001; Wolfram, 2002; Deutsch and Dormann, 2005). To our knowledge, this is the first study in which such a problem is tackled with an antagonistic network involving a large number of interacting strains.

## 2. Materials and Methods

### 2.1. Bacterial Collection

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.

**Figure 1. Antagonism matrix for the bacterial strains considered in this work**. A black box indicates that the bacterial strain in the corresponding column antagonizes that in the box row. Some strains do not appear in the matrix columns because they antagonize no other bacterial strain. The bacterial strains were labeled in decreasing order according to their aggressiveness index AI (numbers of other strains antagonized by it, minus number of strains antagonizing it). In Table 1, we give the original label and the AI number of each strain.

### 2.2. Computational Algorithm

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:

1. 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.

2. The grid is initialized by randomly distributing a fixed number of individual colonies of each bacterial strain among the grid cells.

3. After initializing the grid, the state of all the grid cells is updated according to the following rules:

(a) Let *i* denote the current grid cell. Randomly chose one of its 8 neighboring cells and denote it by *j*.

(b) If both the *i*- and *j*-th grid cells are empty at the *n*-th time step, the *i*-th cell remains empty at time *n* + 1.

(c) If the *i*-th grid cell is empty and the *j*-th grid cell is occupied by a given strain at time *n*, the *i*-th cell is colonized at time *n*+1 with probability *P _{g}*(μ

_{j}) by the strain in the

*j*-th cell, μ

_{j}. Function

*P*(μ

_{g}*) is specified in the Results section for different simulations.*

_{j}(d) If the *i*-th grid cell is occupied at time *n* by strain μ* _{i}*, it becomes empty at time

*n*+ 1 with probability

*P*(μ

_{d}*). The corresponding probability distribution is specified in the Results section for different simulations. If this does not happen, consider the following two possibilities:*

_{i}• If both grid cells (*i* and *j*) are occupied at time *n* and the strain in the *j*-th grid cell antagonizes the strain in the *i*-th grid cell, then the *i*-th grid cell gets empty at time *n* + 1.

• If both grid cells (*i* and *j*) are occupied and no antagonism exists, or if the *i*-th grid cell is occupied and the *j*-th grid cell is empty, the *i*-th cell remains the same.

4. Step 3 is iteratively repeated a fixed number of times.

### 2.3. Generation of Random and Experimental-like Antagonism Matrices

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 construct what we call experimental-like antagonism matrices, we first noted that, according to the original matrix, bacterial strains can be classified as high-level, medium-level and low-level antagonistic strains (a k-means algorithm, MacQueen, 1967 was employed to carry out this classification). Then, we made use of a Markov-Chain Monte Carlo algorithm (Gilks, 2005) to randomly generate *N* links in such a way that the resulting antagonism matrix has the same number of high-level, medium-level, and low-level antagonistic strains as the original matrix.

## 3. Results

### 3.1. Emergency of Community Patches Provides a Survival Mechanism to Low-Level Antagonistic Bacteria

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*(μ) = 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,

_{d}*AI*(μ), equals the number of other strains antagonized by it, minus the number of other strains that antagonize it.

**Figure 2. Emergence of community patches as a result of bacterial growth**. Color indicates strain *Aggressiveness Index* (*AI*). The less aggressive bacteria are represented in light green and the most aggressive ones in maroon; deep blue denotes an empty grid cell. The *AI* of the *i*-th bacterial strain is computed as the number of other strains antagonized by it minus the strain count that antagonize it. **(A)** Initial distribution of the 78 bacterial strains. For each strain 50 initial individual colonies were distributed randomly. **(B)** Final distribution after 500 iteration steps of the simulation algorithm. Total coverage of the surface and emergence of bacterial patches isolating populations can be observed. Beside both graphs we present close-ups of selected regions across the grid, with the ID number of the strain inhabiting each grid cell indicated within (see Table 1).

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.

**Figure 3. Positive long-term correlation between AI and bacterium population. (A)** Box plot summarizing the patch-size statistics of 100 simulations. Bacterial strains are shown in decreasing order according to their AI value. Note that a positive correlation exists between aggressiveness index and average patch size.

**(B)**Correlation between stationary bacterial populations, averaged over 100 simulations, and

*AI*values. A constant was added to the previously defined

*AI*so that it attains positive values for all strains.

**(C)**Growth curves for all the bacterial strains tabulated in Table 1. The AI value corresponding to each strain is indicated by means of the same color code as in Figure 2. Notice that, initially, all strains populations increase monotonically. However, once the grid becomes saturated, antagonistic interactions make the population of low

*AI*strains decrease.

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.

### 3.2. Bacterial-Community Spatial Structure is Disrupted Under Constant Perturbation and all Vulnerable Strains are Driven to Extinction

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).

**Figure 4. Disruption of the emergent community patches by regularly shuffling of the grid. (A)** Initial distribution of the 78 bacterial strains. For each strain 50 initial individual colonies were distributed randomly. **(B)** Final distribution after 500 iteration steps of the simulation algorithm. Every 10 iterations the grid was randomly shuffled. As in the simulations without shuffling, the grid ends up being completely occupied. However, a structure without patches emerges, and all susceptible bacterial strains are missing.

**Figure 5. Susceptible strains become extinct due to regular strain shuffling. (A)** Box plot summarizing the patch-size statistics of 100 simulations. Bacterial strains are shown in decreasing order according to their AI value. Many strains become extinct and the patch size for all the surviving-strain colonies is exactly 1. That is, no patches are formed whatsoever. **(B)** Correlation between *AI* and stationary population levels. Observe that all susceptible strains (most of them with low antagonism levels) become extinct and that the surviving strains occupy the space left by those that become extinct. **(C)** Growth curves for all bacterial strains in this Figure. Notice that regular reshuffling of the grid does not affect bacterial population levels once a stationary state has established.

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.

### 3.3. Influence of the Antagonism Matrix Architecture on the System Dynamic Behavior

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).

### 3.4. Influence of Growth Rate and Initial Population Size on the System Dynamics

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*= 1 and

_{g}*P*= 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).

_{g}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).

### 3.5. Non-Zero Death Rate may Increase the Chance of Survival of Less Antagonistic Bacteria

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*(μ) = 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.

_{d}## 4. Discussion and Conclusions

The computational modeling framework employed in this work has proved to be very useful to studying the spatial and temporal behavior of diverse biological phenomena (Ermentrout and Edelstein-Keshet, 1993; Alber et al., 2003; Rohde, 2005). Hence, we decided to employ it to test the hypothesis that macroscopic community patches can emerge as the result of microscopic individual interaction in a homogeneous environment. Our results not only confirm the feasibility of such hypothesis, but also show that these patches allow bacteria to minimize conflicts while preserving biodiversity.

We tested the robustness of our results by considering different initial-condition scenarios, as well as non-zero death rates and distinct growth rates for different bacterial strains. In all cases we obtained qualitatively equivalent results, thus confirming that the achieved conclusions do not depend on these factors.

On the other hand, it must be emphasized that the rules underlying the implemented computational algorithm are concomitant with low-motility bacteria growing in a homogeneous surface, interacting through pairwise-local-interactions. These restrictions are consistent with the environmental conditions observed in Cuatro Cienegas ponds (Johannesson et al., 2004), where oligotrophic constant conditions have given rise to a great diversity of bacteria who take advantage of any component in the media around them (Escalante et al., 2008; Cerritos et al., 2011) and compete with direct neighbors (Pérez-Gutiérrez et al., 2013; Aguirre-von Wobeser et al., 2014). However, we wondered to what extent having low-motility bacteria is a necessary condition for biodiversity preservation. To test this, we repeated our simulations but included periodical and random shuffling of the grid cells. Since we invariably observed that all susceptible strains become extinct, we conclude that having a slowly changing environment is mandatory for sustaining biodiversity when highly antagonistic, neutral, and highly susceptible strains share the ecosystem.

Previous studies have shown the importance of biodiversity in food-webs, being the web architecture the cause and effect of biodiversity prevalence (Sole and Montoya, 2001; Dunne et al., 2002; Ives and Carpenter, 2007; Allesina and Pascual, 2008), specially under perturbation scenarios (Girvan et al., 2005; Pascual and Guichard, 2005; Dunne and Williams, 2009; Baho et al., 2012). Counter-intuitively our results suggest that antagonism interaction networks may have a similar effect. Previous similar studies have been published (Silvertown et al., 1992; Kerr et al., 2002; Kirkup and Riley, 2004; Károlyi et al., 2005; Walshe, 2006), and some report a large repertoire of possible dynamic behaviors (Silvertown et al., 1992; Károlyi et al., 2005). However, to our knowledge, this is the first study in which a large set of experimental antagonism data is considered.

We are conscious that the model here introduced does not provide a detailed picture of the real-life system. Instead, it can best be regarded as a very simple cartoon or toy model. This is so because the model ignores dynamic aspects that play important ecological roles. For instance, positive interactions are well-documented in the case of biofilms, which tend to aggregate various types of bacteria and promote positive interactions among them (in this way, positive interactions can generate micro-habitats which introduce a level of physicochemical heterogeneity even in an otherwise rather stable and constant habitat). In this respect, we build our model following Einstein's advice that every theory (model) should be as simple as possible, but not simpler (i.e., not so simple that it does not represent reality any longer). To our consideration, given the amount of available experimental information, the present is the simplest possible model one can come out with to tackle the question of whether the antagonism matrix found by Pérez-Gutiérrez et al. (2013) can explain a heterogeneous bacterial community distribution in a homogeneous habitat.

## Author Contributions

GO and MS conceived the project, designed the *in-silico* experiments, and interpreted data and results. RZ conducted the experiments and analyzed and interpreted data. All authors contributed to writing the manuscript.

## Conflict of Interest Statement

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.

## Acknowledgments

The authors are grateful to Profs. Michael Travisano, Valeria Souza, Michael Mackey, and Frédéric Guichard for fruitful and encouraging discussions.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fmicb.2015.00489/abstract

## References

Abriouel, H., Franz, C. M., Omar, N. B., and Gálvez, A. (2011). Diversity and applications of bacillus bacteriocins. *FEMS Microbiol. Rev*. 35, 201–232. doi: 10.1111/j.1574-6976.2010.00244.x

Aguirre-von Wobeser, E., Soberón-Chávez, G., Eguiarte, L. E., Ponce-Soto, G. Y., Vázquez-Rosas-Landa, M., and Souza, V. (2014). Two-role model of an interaction network of free-living γ-proteobacteria from an oligotrophic environment. *Environ. Microbiol*. 16, 1366–1377. doi: 10.1111/1462-2920.12305

Alber, M. S., Kiskowski, M. A., Glazier, J. A., and Jiang, Y. (2003). “On cellular automaton approaches to modeling biological cells,” in *Mathematical Systems Theory in Biology, Communications, Computation, and Finance*, eds J. Rosenthal and D. S. Gilliam (New York, NY: Springer), 1–39. doi: 10.1007/978-0-387-21696-6_1

Allesina, S., and Pascual, M. (2008). Network structure, predator–prey modules, and stability in large food webs. *Theor. Ecol*. 1, 55–64. doi: 10.1007/s12080-007-0007-8

Baho, D. L., Peter, H., and Tranvik, L. J. (2012). Resistance and resilience of microbial communities–temporal and spatial insurance against perturbations. *Environ. Microbiol*. 14, 2283–2292. doi: 10.1111/j.1462-2920.2012.02754.x

Bohannan, B. J., Kerr, B., Jessup, C. M., Hughes, J. B., and Sandvik, G. (2002). Trade-offs and coexistence in microbial microcosms. *Antonie van Leeuwenhoek* 81, 107–115. doi: 10.1023/A:1020585711378

Cadotte, M. W. (2007). Competition-colonization trade-offs and disturbance effects at multiple scales. *Ecology* 88, 823–829. doi: 10.1890/06-1117

Cerritos, R., Eguiarte, L. E., Avitia, M., Siefert, J., Travisano, M., Rodríguez-Verdugo, A., et al. (2011). Diversity of culturable thermo-resistant aquatic bacteria along an environmental gradient in cuatro ciénegas, coahuila, méxico. *Antonie Van Leeuwenhoek* 99, 303–318. doi: 10.1007/s10482-010-9490-9

Costello, E. K., Stagaman, K., Dethlefsen, L., Bohannan, B. J., and Relman, D. A. (2012). The application of ecological theory toward an understanding of the human microbiome. *Science* 336, 1255–1262. doi: 10.1126/science.1224203

Czárán, T. L., Hoekstra, R. F., and Pagie, L. (2002). Chemical warfare between microbes promotes biodiversity. *Proc. Natl. Acad. Sci*. 99, 786–790. doi: 10.1073/pnas.012399899

Deutsch, A., and Dormann, S. (2005). Cellular automaton modeling of biological pattern formation. *FASEB* 23, 12. doi: 10.1007/s10710-006-9021-7

Dunne, J. A., and Williams, R. J. (2009). Cascading extinctions and community collapse in model food webs. *Philos. Trans. R. Soc. B Biol. Sci*. 364, 1711–1723. doi: 10.1098/rstb.2008.0219

Dunne, J. A., Williams, R. J., and Martinez, N. D. (2002). Network structure and biodiversity loss in food webs: robustness increases with connectance. *Ecol. Lett*. 5, 558–567. doi: 10.1046/j.1461-0248.2002.00354.x

Dunson, W. A., and Travis, J. (1991). The role of abiotic factors in community organization. *Am. Nat*. 138, 1067–1091. doi: 10.1086/285270

Eisenhauer, N., Scheu, S., and Jousset, A. (2012). Bacterial diversity stabilizes community productivity. *PLoS ONE* 7:e34517. doi: 10.1371/journal.pone.0034517

Ermentrout, G. B., and Edelstein-Keshet, L. (1993). Cellular automata approaches to biological modeling. *J. Theor. Biol*. 160, 97–133. doi: 10.1006/jtbi.1993.1007

Escalante, A. E., Eguiarte, L. E., Espinosa-Asuar, L., Forney, L. J., Noguez, A. M., and Souza Saldivar, V. (2008). Diversity of aquatic prokaryotic communities in the cuatro cienegas basin. *FEMS Microbiol. Ecol*. 65, 50–60. doi: 10.1111/j.1574-6941.2008.00496.x

Evans, M. R., Bithell, M., Cornell, S. J., Dall, S. R., Díaz, S., Emmott, S., et al. (2013). Predictive systems ecology. *Proc. R. Soc. B Biol. Sci*. 280:20131452. doi: 10.1098/rspb.2013.1452

Foster, K. R., and Bell, T. (2012). Competition, not cooperation, dominates interactions among culturable microbial species. *Curr. Biol*. 22, 1845–1850. doi: 10.1016/j.cub.2012.08.005

Gardner, M. (1970). Mathematical games: the fantastic combinations of john conway new solitaire game. *Sci. Am*. 223, 120–123. doi: 10.1038/scientificamerican1070-120

Girvan, M., Campbell, C., Killham, K., Prosser, J., and Glover, L. (2005). Bacterial diversity promotes community stability and functional resilience after perturbation. *Environ. Microbiol*. 7, 301–313. doi: 10.1111/j.1462-2920.2005.00695.x

Hogeweg, P. (1988). Cellular automata as a paradigm for ecological modeling. *Appl. Math. Comput*. 27, 81–100. doi: 10.1016/0096-3003(88)90100-2

Hooper, D., Chapin, F. III., Ewel, J., Hector, A., Inchausti, P., Lavorel, S., et al. (2005). Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. *Ecol. Monogr*. 75, 3–35. doi: 10.1890/04-0922

Ives, A. R., and Carpenter, S. R. (2007). Stability and diversity of ecosystems. *Science* 317, 58–62. doi: 10.1126/science.1133258

Iwasa, Y., Nakamaru, M., and Levin, S. A. (1998). Allelopathy of bacteria in a lattice population: competition between colicin-sensitive and colicin-producing strains. *Evol. Ecol*. 12, 785–802. doi: 10.1023/A:1006590431483

Jessup, C. M., Kassen, R., Forde, S. E., Kerr, B., Buckling, A., Rainey, P. B., et al. (2004). Big questions, small worlds: microbial model systems in ecology. *Trends Ecol. Evol*. 19, 189–197. doi: 10.1016/j.tree.2004.01.008

Johannesson, K. H., Cortés, A., and Kilroy, K. C. (2004). Reconnaissance isotopic and hydrochemical study of cuatro ciénegas groundwater, coahuila, méxico. *J. South Am. Earth Sci*. 17, 171–180. doi: 10.1016/j.jsames.2004.01.002

Károlyi, G., Neufeld, Z., and Scheuring, I. (2005). Rock-scissors-paper game in a chaotic flow: the effect of dispersion on the cyclic competition of microorganisms. *J. Theor. Biol*. 236, 12–20. doi: 10.1016/j.jtbi.2005.02.012

Kerr, B. (2007). “The ecological and evolutionary dynamics of model bacteriocin communities,” in *Bacteriocins* eds A. R. Margaret and A. C Milind (Berlin; Heidelberg: Springer), 111–134. doi: 10.1007/978-3-540-36604-1_6

Kerr, B., Riley, M. A., Feldman, M. W., and Bohannan, B. J. (2002). Local dispersal promotes biodiversity in a real-life game of rock–paper–scissors. *Nature* 418, 171–174. doi: 10.1038/nature00823

Kirkup, B. C., and Riley, M. A. (2004). Antibiotic-mediated antagonism leads to a bacterial game of rock–paper–scissors *in vivo*. *Nature* 428, 412–414. doi: 10.1038/nature02429

Kneitel, J. M., and Chase, J. M. (2004). Trade-offs in community ecology: linking spatial scales and species coexistence. *Ecol. Lett*. 7, 69–80. doi: 10.1046/j.1461-0248.2003.00551.x

Kuznetsov, S., Dubinina, G., and Lapteva, N. (1979). Biology of oligotrophic bacteria. *Annu. Rev. Microbiol*. 33, 377–387. doi: 10.1146/annurev.mi.33.100179.002113

MacQueen, J. (1967). “Some methods for classification and analysis of multivariate observations,” in *Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability*, Vol. 1 (Berkeley, CA), 281–297.

Momeni, B., Chen, C.-C., Hillesland, K. L., Waite, A., and Shou, W. (2011). Using artificial systems to explore the ecology and evolution of symbioses. *Cell. Mol. Life Sci*. 68, 1353–1368. doi: 10.1007/s00018-011-0649-y

Pace, M. L. (1988). Bacterial mortality and the fate of bacterial production. *Hydrobiologia* 159, 41–49. doi: 10.1007/BF00007366

Pascual, M., and Guichard, F. (2005). Criticality and disturbance in spatial ecological systems. *Trends Ecol. Evol*. 20, 88–95. doi: 10.1016/j.tree.2004.11.012

Pérez-Gutiérrez, R.-A., López-Ramírez, V., Islas, Á., Alcaraz, L. D., Hernández-González, I., Olivera, B. C. L., et al. (2013). Antagonism influences assembly of a bacillus guild in a local community and is depicted as a food-chain network. *ISME J*. 7, 487–497. doi: 10.1038/ismej.2012.119

Prasad, S., Manasa, P., Buddhi, S., Singh, S. M., and Shivaji, S. (2011). Antagonistic interaction networks among bacteria from a cold soil environment. *FEMS Microbiol. Ecol*. 78, 376–385. doi: 10.1111/j.1574-6941.2011.01171.x

Prosser, J. I., Bohannan, B. J., Curtis, T. P., Ellis, R. J., Firestone, M. K., Freckleton, R. P., et al. (2007). The role of ecological theory in microbial ecology. *Nat. Rev. Microbiol*. 5, 384–392. doi: 10.1038/nrmicro1643

Raes, J., and Bork, P. (2008). Molecular eco-systems biology: towards an understanding of community function. *Nat. Rev. Microbiol*. 6, 693–699. doi: 10.1038/nrmicro1935

Rohde, K. (2005). Cellular automata and ecology. *Oikos* 110, 203–207. doi: 10.1111/j.0030-1299.2005.13965.x

Sarkar, P. (2000). A brief history of cellular automata. *ACM Comput. Surv*. 32, 80–107. doi: 10.1145/349194.349202

Servais, P., Billen, G., and Rego, J. V. (1985). Rate of bacterial mortality in aquatic environments. *Appl. Environ. Microbiol*. 49, 1448–1454.

Silvertown, J., Holtier, S., Johnson, J., and Dale, P. (1992). Cellular automaton models of interspecific competition for space–the effect of pattern on process. *J. Ecol*. 80, 527–533. doi: 10.2307/2260696

Sole, R. V., and Montoya, M. (2001). Complexity and fragility in ecological networks. *Proc. R. Soc. Lond. Ser. B Biol. Sci*. 268, 2039–2045. doi: 10.1098/rspb.2001.1767

Souza, V., Espinosa-Asuar, L., Escalante, A. E., Eguiarte, L. E., Farmer, J., Forney, L., et al. (2006). An endangered oasis of aquatic microbial biodiversity in the chihuahuan desert. *Proc. Natl. Acad. Sci. U.S.A*. 103, 6565–6570. doi: 10.1073/pnas.0601434103

Walshe, R. (2006). “Modeling bacterial growth patterns in the presence of antibiotic,” in *Engineering of Complex Computer Systems, 2006. ICECCS 2006. 11th IEEE International Conference on* (Piscataway, NJ: IEEE), 10.

Keywords: bacterial antagonism, ecological modeling, community emergence, spatial patterns, bacterial biodiversity, cuatro cienegas

Citation: Zapién-Campos R, Olmedo-Álvarez G and Santillán M (2015) Antagonistic interactions are sufficient to explain self-assemblage of bacterial communities in a homogeneous environment: a computational modeling approach. *Front. Microbiol*. **6**:489. doi: 10.3389/fmicb.2015.00489

Received: 16 October 2014; Accepted: 04 May 2015;

Published: 21 May 2015.

Edited by:

Michael Travisano, University of Minnesota, USAReviewed by:

Sébastien Duperron, Université Pierre et Marie Curie, FranceElisabeth Margaretha Bik, Stanford University School of Medicine, USA

Copyright © 2015 Zapién-Campos, Olmedo-Álvarez and Santillán. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Moisés Santillán, Unidad Monterrey, Centro de Investigación y Estudios Avanzados del IPN, Vía del Conocimiento 201, Parque PIIT, Apodaca, 66600 México, msantillan@cinvestav.mx