Modeling evolution of spatially distributed bacterial communities: a simulation with the haploid evolutionary constructor

Background Multiscale approaches for integrating submodels of various levels of biological organization into a single model became the major tool of systems biology. In this paper, we have constructed and simulated a set of multiscale models of spatially distributed microbial communities and study an influence of unevenly distributed environmental factors on the genetic diversity and evolution of the community members. Results Haploid Evolutionary Constructor software http://evol-constructor.bionet.nsc.ru/ was expanded by adding the tool for the spatial modeling of a microbial community (1D, 2D and 3D versions). A set of the models of spatially distributed communities was built to demonstrate that the spatial distribution of cells affects both intensity of selection and evolution rate. Conclusion In spatially heterogeneous communities, the change in the direction of the environmental flow might be reflected in local irregular population dynamics, while the genetic structure of populations (frequencies of the alleles) remains stable. Furthermore, in spatially heterogeneous communities, the chemotaxis might dramatically affect the evolution of community members.


Background
Prokaryotes considered as the most ancient living organisms and essential part of the Earth biosphere. In particular, prokaryotic (or microbial) communities maintain all major biogeochemical cycles [1][2][3]. Typical examples of the communities are spatially complex, layered structures of the bacterial mats or biofilms [4][5][6]). A majority of prokaryote species cannot be cultured, so we have to study them within their natural environment, i.e. in communities. Hence, mathematical modeling and simulation of bacterial communities are indispensable for understanding of the functioning and evolution of bacteria. Spatial factors are well-known to be among major forces of the evolution [7][8][9][10][11]. From ecological and evolutionary points of view, spatial distribution of species plays a large role in local microbial cooperation and competition; consequently, spatial distribution influences evolution [4,[12][13][14][15]. Combined with other evolutionary factors, these factors affect dynamics of allele frequencies in populations of a community [16][17][18][19][20]. Thus, it has been shown that mutator populations adapted faster than wild-type populations in both liquid and solid environments. Also, it has been shown that independently of the mutation rate, the increase in fitness in the spatially structured environment was smaller than in the unstructured one [21]. Summarizing the mentioned above, the study of dependence of functioning and evolution of microbial communities on the spatial distribution of organisms and substances is of essential theoretical interest. Equally important is the study of how various communities transform their habitats.
At present, a number of software tools are available for modeling and simulation of spatially distributed microbial communities. A majority of these tools, such as cellular automata UMCCA [22], multi-agent software packages AgentCell [23], AQUASIM [24], INDISIM [25] and others [26][27][28][29][30], emphasize the details of the spatial distribution of cells as such. However, the study of the genetic variability effects upon a spatial structure of a community is of equal importance. On the other hand, evolutionary oriented software packages, for example, AEvol [31,32] and others [33][34][35], mainly focus on genetic structure, evolution and/or metabolism and do not provide the modules for the study of spatial organization. In present study, we took into account both the spatial distribution of cells and substances and the genetic variability of individual populations by developing a software with expanded set of options that permits the requested modeling.
In this paper, we have carried out the modeling of the spatial heterogeneity of environment surrounding the microbial community and its influence on population dynamics and evolution. For this purpose, we have developed a software package for the modeling of spatially distributed microbial communities (1D, 2D, and 3D versions); the package was added to Haploid Evolutionary Constructor (HEC) software that was described previously [36]. We have constructed and analyzed computational models of the prokaryotic communities' evolution with "poisoner-prey" trophic interactions in heterogeneous flow environments. We have also studied models of the horizontal gene transfer in prokaryotic communities living in spatially distributed habitats under changing environmental conditions. According to our models, a combination of chemotaxis with other spatial factors might significantly affect the life of prokaryotic communities through the changes in both the dynamics of population and its genetic structure.

Modeling methods
Simulations have been carried out with the HEC software. The key object modeled in the HEC is the prokaryotic polymorphic population which is assumed here to be equivalent to species (or strain). Populations consist of cells utilizing substrates which are either consumed from the environment or are synthesized by themselves. Utilization energy is then used for the reproduction and synthesis of other substrates. Synthesized substrates may be either used by cells for own requirements or secreted into habitat. In the latter case, substrates may be consumed by cells of other species. In the HEC, substrate synthesis and utilization, as well as cell reproduction, are described via the corresponding gene networks (GN) which are named strategies (i.e. there are synthesis strategies and reproduction strategies, see details in Figure 1). Numerical parameters of those GNs assumed to be "genes" inheriting across the generations. Cells are supposed to belong to the same population (species) if they possess the same GN structures. However, numerical parameters of GNs in cells of a certain population may vary that models the genetic polymorphism. The generalized genome of a polymorphic population is formed by a set of allele distributions for each "gene". Mutations changing numerical values of GN parameters consequently change the corresponding allele distribution. Mutations may either be pre-described by the user, or be randomly generated with beta distribution (the user may also control its parameters). In the HEC, both finite and infinite sites models may be used optionally. At the same time, horizontal gene transfer or genes changes the structure of GN, which, in turn, aids in generation of novel species. These processes may also be pre-described by the user or be randomly generated. Thus, the polymorphic population is characterized by the generalized population genome, the synthesis strategy, the reproduction strategy, intercellular substrates and other parameters (Additional file 1 fig. S1, detailed description is published in [37,38]). Most of the HEC parameters concerning the cell size, amounts of substrates required for cell reproduction and other factors have been estimated on the base of E.coli data [39]. The simulations describe the dynamics of populations, dynamics of allele frequencies (due to either selection or mutations), origin and extinction of species, and dynamics of environmental substrates. All the dynamics of the system is conditioned by the efficiency of metabolic processes of community members in particular environmental conditions (which in turn depend on corresponding GN). Number of classical models including the logistic growth [40] of the bacterial population and Fisher's fundamental theorem of natural selection [41,42] had previously been implemented in the HEC and found to be consistent with the model expectations [36]. We have extended the original HEC (uniform mixing case, denoted as 0D [36]) by adding the cases of one-, two-or three-dimensional spatial distributions of cells and substrates into the model (1D, 2D and 3D versions of the HEC, correspondingly). We use the grid of "point environments" -finite mesh of little volume connected with continuous flow (Figure 2) to model spatial distributions in HEC 1D-3D.
Standard simulation step in the HEC 1D-3D consists of the following two stages: (1) calculation of new states for each point environments (which is independent and can be performed simultaneously) including simulation of the following processes: consumption of substrates, utilization of substrates, reproduction, substrate synthesis and secretion. This stage is apparently inherited from the HEC 0D; (2) spatial redistribution of substrates and cells in the modeled system including the simulation of flow, diffusion, and chemotaxis processes ( Figure 2). Detailed description of spatial processes, as well as nuances related to differences in characteristic times of reproduction and spatial processes, is presented in Additional file 1 ( fig. S3, S4).
Hence, the processes of substrate production and utilization, reproduction, mutation, genes loss and horizontal transfer are simulated as part of the standard HEC 0D iteration and occur independently in the each mesh point (Additional file 1 fig. S2). Only spatial redistribution of organisms and substrates requires nodes synchronization (Additional file 1 fig. S3).

Results and discussion
Modeling "poisoner-prey" community in spatially heterogeneous habitats Using the methods described above, we have studied the change of genetic diversity in populations of the "poisoner-prey" community. The "poisoner-prey" model has been previously described in our recent studies [36,43]. One-dimensional flow-through environment (1D tube, Figure 3). Microbial community consisted of two populations -the poisoner and the prey ( Figure 4). Metabolic by-product secreted by the poisoners would inhibit the growth of the prey which product would conversely activate the growth of the poisoner population. Besides, both populations would consume non-specific substrate which had come into the habitat with the inflow. They used it for the growth. In this model (as well as in the others in this paper) we have considered systems of various nodes  number (from 10 to 1000). Since a qualitative character of model behavior remained the same, we have considered 10-node tubes without losing generality in this paper. Furthermore, as we have not simulated any mutation, the allele frequencies changed only by the efficiency of metabolic processes of community members in particular environment, which, in turn, depends on corresponding GN). For each allelic combination in a polymorphic population P, the fitness is calculated as the (P after -P before )/P before ratio, where P before and P after denote the size of P before and after the reproduction process, respectively.
In the described model, we have considered the following ways of a spatial organization: First, we have considered a model (Additional file 2) with the uniform initial distribution of cells and substances, and flow-through ( Figure 3a) spatial organization with the flow rate of 0.02 (i.e. 2% of cells and substances of the node is taken away at each iteration). Initial populations were genetically polymorphic: prey cells varied in terms of their sensitivity to the S2 inhibitor, poisoner cells varied in terms of their efficiency of the S1 substrate utilization. Population dynamics of the community demonstrates (Additional file 1 Fig. S6) that after a relatively short period of oscillation in the preys' size which is associated with the steady growth of the poisoners, the size of both populations becomes stable.
However, if we look at separate nodes ( Figure 5), we can see that after 250 generations the prey survives only in 1 st and 2 nd nodes, while the poisoner lives in all nodes of the habitat. We think that this is due to the preys which live only on the non-specific substrate which is in sufficient amount only in the nodes close to the inflow. At the same time, the poisoner can partially compensate the lack of the non-specific substrate with the S1.
At the same time, we have analyzed the change of genetic diversity in both populations. Figure 6 shows   that adaptive alleles (i.e. in case with the prey, these are parameters of sensitivity to S2, in case with the poisoner, these are efficiencies of S1 utilization) completely displace less adaptive ones in course of time. Moreover, the nodes outlying from the source of the non-specific substrate are characterized by the rapid genetic diversity loss for the prey population which is associated with the preservation of genetic diversity in the poisoner (Table 1 node 10 in Figure 6). Reverse situation is observed in proximal nodes (Table 1 node 1 in Figure 6).
However, in central nodes (node 5 in Table 1) we have observed an increase of time of the genetic diversity preservation -inadaptive alleles would extinct only at 104 ± 0.3 generation. We think this can be explained by the flow transferring migrants from proximal regions of a habitat in which relatively high genetic diversity remains.
Therefore, spatial localization of microorganisms may influence the evolutionary rate: depending on the cell position relative to the source of non-specific substrate, we have observed evolutionary rates among poisoners and preys to be differed from each other. Poisoners evolve more rapidly than preys when located near the source of non-specific substrate and vice versa. Located far from the source, the preys evolve more rapidly.
Later we have studied the above model by adding a nonuniform initial state and chemotaxis. Spatial structure of a habitat has been set to perpendicular-flow (Figure 3b). Initial distribution of the inhibitor S2 has been set to gradient (the highest concentration was in the node 1, the lowest in the node 10). Flow rate was the same as in the previous case. We have analyzed the influence of chemotaxis on the model behavior. Figure 7 shows the chemotaxis-off case (Additional file 3). In different nodes oscillations in the prey population size vary due to nonuniform distribution of inhibitor. Poisoners population size negligibly varies in different nodes. Situation changes when chemotaxis switches on: cells can actively move to nodes with more favorable conditions. Simulation results are shown in Additional file 1 fig. S7, S8. Inhibitor gradient and cell movement cause an origin of local oscillations in population dynamics which affects the dynamics of both separate nodes and a habitat in whole.  Node 10 582 ± 2.9 Node 10 19 ± 0 One hundred simulation runs have been performed in which the initial allele distributions varied in such a way that the frequency for each allele lied in the range 5-15% (total number of alleles in a distribution was 9). Confidence intervals are stated at the 95% confidence level. Analysis of this model under various environmental conditions has shown that the non-uniform distribution of substrates that was due to either initial inequalities data or high diffusion rate along with the chemotaxis resulted in emergent community dynamics evident in a movement of cells to a better place of living, into the nodes with higher concentrations of necessary substrates. Under these circumstances, both local oscillations and even an irregular behavior were observed (Additional file 1 fig. S9). This finding is in an agreement with experimentally observed behavior of antagonistic E.coli strains reported earlier [44].
Detailed analysis of the genetic structure dynamics in these populations is described below. We have considered the "poisoner-prey" model with a perpendicular flow and a diffusion-dependent non-uniformity of substrates. At the baseline, the genetic polymorphism was same to that in previous models ( Figure 6). Population dynamics of the model is shown in Additional file 1 fig. S10.
Diffusion causes non-uniformity. Combined with active movement of cells, it leads to irregularities in the dynamics within local populations (Additional file 1 fig. S10). It is evident that chemotaxis may promote irregular oscillations of the population size in separate nodes while at the genetic level the dynamics of allele frequencies remains stable (Additional file 1 fig. S11).

Modeling horizontal transfer of genes in spatially distributed systems under varying environmental conditions
The second model describes co-functioning and competition of populations in a community consisting of two trophic cycles living in the 1D tube habitat (Figure 3) under varying environmental conditions. This model is the extension of a previously published model with an uniform mixing (0D) [38]. Each trophic cycle includes three populations (Figure 8). For this community, we have considered two combinations of breeding strategies: NC-NC and C-C. Here, NC-NC strategy implies that all populations employ non-compensatory strategy, implementing the Liebig's law of the minimum: "The growth is controlled not by the total amount of resources available, but by the most scarce resource (limiting factor)" [45]; C-C means that all populations of a community breed via the compensatory reproduction strategy, implementing the Rubel's law of compensation of ecological factors: "If one factor intensifies the action of another, the minimum of the latter is less than it would be without the help of the former. The contrary effect is possible when the accompanying factor results in raising the minimum" [45] (formulas for both strategies are presented in the Supplementary materials). Our previous studies for the 0D case have shown that during long period of starvation the non-compensatory strategy preserved the genetic diversity while compensatory strategy was preserving it only for a shortwhile [38]. At the same time, non-compensatory systems have been shown to be less adaptable for harsh environmental conditions than compensatory ones.
We have considered both cases with and without chemotaxis. For these models, we have studied how the horizontal transfer (HT) of genes and change of environmental conditions affected the community functioning: 1. On the 100 th generation, we have simulated the HT of the gene of the S6 specific substrate utilization from cells of the P6 population into cells of the P1 population. It would lead to the origin of new type cells which then would form the P7 population connecting two trophic cycles. 2. On the 2000 th generation, we have simulated the decrease of a non-specific substrate concentration in the flow by factor of 100. Such a "starvation" continued for 1000 generations and then initial conditions had been restored.
First, we have considered flow-through habitats. In C-C communities when chemotaxis was switched off, the HT would change the fate of the communities if only occurred in nodes close to the source of a nonspecific substrate (nodes 1-4, Additional file 1 fig. S12): new population would survive and population sizes would differentiate. If the HT occurred in outlying nodes (5)(6)(7)(8)(9)(10), new population could not settle and would be eliminated soon. Chemotaxis encouraged fixation (Additional file 4 and Additional file 5) and sustaining of biodiversity (surviving of the P7 population). However, total biomass of such a community turned out to be lower (Figure 9). Palaeontologist Zherikhin and paleobotanist Meyen proposed the concepts of zonal stratification [46] and phytospreading [47], which are the extensions of Darlington's "equatorial pump" concept [48]. According to these concepts, the most rates of taxa formation are observed in such locations in which the struggle for existence against abiotic factors is reduced. In the previous model (chemotaxis off case), the P7 population survived only if the HT occurred in substrate-rich nodes (i.e. close to the node 1), which are analogues of Zherikhin's and Meyen's biotopes. In poor ecosystems, the novel species could not grow up to the necessary size for an effective competition in the community in spite of the fact that novel genotype would be potentially more adaptive than others (more substrates could be utilized). It is in good agreement with the concepts of equatorial pump and phytospreading, as the evolutionary success of a novel species is related to not only preadaptations, but also to habitat conditions [49]. There is also an interesting conclusion of the habitat constraint mentioned above. It is commonly supposed that adaptations appear in one or several individuals. If they occur in unfavorable conditions, a small subpopulation of mutants would not fixate in a community owing to its size. That is why pre-adaptations should be accumulated in the form of a neutral variability during the period of relatively optimal habitat conditions. They manifest when conditions change.
Simulation results were similar to those for compensatory communities in non-compensatory communities without chemotaxis. However, in contrast to the C-C case, the presence of chemotaxis has dramatically changed the fate of NC-NC communities. HT has destabilized the system leading to the extinction of the part of the acceptor trophic ring (Figure 10): if the HT occurred close to the source of a non-specific substrate (Additional file 6), the whole acceptor trophic ring would perish ( Figure 11); in other cases (Additional file 7 and Additional file 8) only the P2 population would perish ( Figure 12). Finally, we have considered systems with perpendicularflow habitats. In C-C communities, chemotaxis was found to result in a bit more expressed differentiation of the population size under starvation (Additional file 1 Figure 8 Trophic graph of a community consisting of two trophic cycles: P1-P2-P3 and P4-P5-P6 [38]. N is a non-specific substrate contained in the flow. HT of the gene from P6 cells to P1 cells is shown by a lightning. As a result of the HT, the new type of cells, which forms P7 population (grey arrow) originates.
fig. S13) -P max /P min ratios are 1.77 (chemotaxis is off, Additional file 9 ) and 1.97 (chemotaxis is on, Additional file 10). In NC-NC communities (Additional file 11 and Additional file 12 ), the cell movement ability would lead to the origin of irregular oscillations (Additional file 1 fig. S14), although of a severely limited amplitude (~4*10 7 cells).
It seems surprising to us that while the P2 is the potential beneficiary of the community perturbations caused by the HT (initially the P2 cells were fed only by P1 cells, and after HT -by both P1 and P7). It should particularly be noted that the destabilization of the community occurred prior to the starvation. The novel P7 population is rather adaptive as it always survives. However, the P7 does not dominate in the community, which is rather not obvious. In our opinion, the reason is that the NC trophic ring consists of the highly specialized symbiotrophic populations realizing the optimal system of trophic interactions. If the HT along with the chemotaxis promotes the origin of high-competing species and then such species may destabilize the established community even long before the starvation. It consequently may lead to unpredictable dynamics of the community.
The location of the HT significantly affects the fate of communities living in flow-through habitats leading to the differentiation and even the change of the community structure. In perpendicular-flow habitats, the location of the HT plays much lesser role: functional modes of a community do not change whenever the HT occurred.

Conclusion
In this study, the methods for modeling spatially distributed microbial communities with changing genetic structure have been presented along with the HEC 1D-3D software package. The software allows one to build a model of the microbial community that takes into  account both spatial environmental factors and genetic and metabolic variations. Thus, the software described provides it users with expanded simulation capabilities and is superior to existing tools for simulation of microbial communities.
The models of functioning and evolution of prokaryotic communities implementing "poisoner-prey" trophic interactions in spatially heterogeneous habitats have been constructed and analyzed. Spatial localization of organisms has been shown to affect the selection intensity and their evolutionary rate: depending on the distance from a substrate source, we have observed evolutionary rates among poisoner and prey populations to be differed from each other. For "poisoner-prey" communities living in perpendicular-flow habitats, we have also shown that spatial heterogeneity might lead to the origin of irregular population dynamics, but dynamics of genetic structure (allele frequencies) at the time would remain stable.
The horizontal transfer models in symbiotic communities consisted of two trophic cycles and living in spatially distributed habitats under varying environmental conditions have also been analyzed. Location of the HT origin has been noted to essentially influence the fate of the community, especially in flow-through habitats. As for the temporary starvation, various locations of the HT have determined various population structure and even species composition of the community. In perpendicular-flow habitats, place of the HT origin would play much lesser role: functional regimes of the community would not change whenever it occurred. Therefore, it has been shown that movement factors (chemotaxis) associated with a non-uniform spatial distribution of cells and substances might dramatically affect the life of communities which has been manifested in both population dynamics and dynamics of genetic diversity of populations for all considered models.