Agent-based modeling of competence phenotype switching in Bacillus subtilis
© Stiegelmeyer and Giddings; licensee BioMed Central Ltd. 2013
Received: 9 November 2012
Accepted: 21 March 2013
Published: 3 April 2013
It is a fascinating phenomenon that in genetically identical bacteria populations of Bacillus subtilis, a distinct DNA uptake phenotype called the competence phenotype may emerge in 10–20% of the population. Many aspects of the phenomenon are believed to be due to the variable expression of critical genes: a stochastic occurrence termed “noise” which has made the phenomenon difficult to examine directly by lab experimentation.
To capture and model noise in this system and further understand the emergence of competence both at the intracellular and culture levels in B. subtilis, we developed a novel multi-scale, agent-based model. At the intracellular level, our model recreates the regulatory network involved in the competence phenotype. At the culture level, we simulated growth conditions, with our multi-scale model providing feedback between the two levels.
Our model predicted three potential sources of genetic “noise”. First, the random spatial arrangement of molecules may influence the manifestation of the competence phenotype. In addition, the evidence suggests that there may be a type of epigenetic heritability to the emergence of competence, influenced by the molecular concentrations of key competence molecules inherited through cell division. Finally, the emergence of competence during the stationary phase may in part be due to the dilution effect of cell division upon protein concentrations.
The competence phenotype was easily translated into an agent-based model – one with the ability to illuminate complex cell behavior. Models such as the one described in this paper can simulate cell behavior that is otherwise unobservable in vivo, highlighting their potential usefulness as research tools.
KeywordsAgent-based modeling Bacillus subtilis Bistable switch Competence
The competence state exhibited by Bacillus subtilis is an example of a bacterial phenotype driven by changes in gene and protein expression states rather than changes in genotype. Competence is a DNA uptake mechanism that appears to be a cell survival strategy for either procuring new genetic information or obtaining DNA as food. The emergence of competence is correlated with high cell density and nutrient limiting conditions. In those conditions, approximately 10–20% of a B. subtilis population will express the competence phenotype. The mechanism by which a small fraction of the population becomes competent is presently attributed to the variable expression of specific genes, i.e., genetic “noise.” “Noise” in gene expression at present is poorly understood.
Stochastic intracellular molecular interactions and environmental inputs like changing nutrients and cell density are known to initiate phenotype switching and are sources of noise[3, 10–12]. Yet, the precise molecular interactions driving variable gene expression are difficult to capture in a typical laboratory environment. To educate ourselves about the molecular interactions involved in the competence phenotype, and to speculate if simple temporal and spatial molecular interactions also contribute to genetic noise, we developed a virtual model of a simplified B. subtilis cell in a cell colony-like environment.
A variety of mathematical modeling methods have been employed in an attempt to understand the nature of the competence switching process[2, 13–15]. These models have addressed the stochastic nature of competence by modeling noise in the system with varying degrees of specificity, using pre-defined noise terms and the Gillespie stochastic modeling algorithm. However, to model noise, a stochastic process, we felt that employing a modeling technique that is itself inherently stochastic may offer additional insights about the nature of competence.
We decided to employ an alternative modeling technique – agent-based modeling (ABM). We hypothesized that the randomly interacting agents of an ABM would be suitable for modeling the genetic noise observed in the competence switch. In our ABM, agents simulate large molecules like proteins, RNAs, and DNAs. Agents interact with each other and their environment based on a set of well-defined rules, either literature-based or biologically intuitive.
Recently, ABMs have been used with great success to model phenomena in bacteriology such as biofilm development[17–21], the transmission dynamics of antibiotic resistance, and antibiotic resistance mechanisms in Staphylococcus aureus. We recently reported the use of ABMs as applied to the chemotactic switching system in Escherichia coli.
Here we report the ABM we designed to mimic the B. subtilis competence phenotype. The model was tuned to simulate a colony of B. subtilis cells where competence emerged in 10–20% of the population. We make the following observations and predictions about genetic “noise”: a) random spatial-temporal interactions lead to competence, b) daughter cells which inherit portions of the competence machinery are more likely to exhibit the competence phenotype, and c) dilution events like cell division inhibit competence emergence until stationary phase. The resulting model’s source code is openly available for exploration at http://www.giddingslab.org/software, for use with the open-source modeling platform Repast Simphony version 1.2.
Modeling environment and overview
Agents, the rules they execute, and the agents they interact with
ComK, ComX, DegU, Repressor
move, bind, translation
Promoter, ComK, MecA
ComK, ComS, ClpC/ClpP
move, bind, death
Nutrients, ComX Peptides
move, shove, generatePeptide, consumePeptide, consumeNutrients, life, death
Nutrients, ComX Peptides
Time in the model was represented as discrete updates to the state of all agents in the model. A time step was completed when all agent rules had been attempted in a random order. Rules for both the cellular and intracellular models are handled in this manner, keeping both models on the same time scale.
Nutrients and the comX peptide were mathematically modeled using diffusion equations adapted to the ABM environment (Figure 2). The model uses continuous equations due to the high concentrations of nutrients and eventual comX. Concentrations are monitored in each cell of the grid environment and diffusion to adjacent cells is calculated by the equations described below.
In the model, agent movement and interactions with other agents are defined by rules. Rule execution is stochastic and subject to meeting a probability threshold after a random draw from the uniform distribution (see Parameter Estimation).
For example, if two molecules have been shown to bind biologically with a high affinity, then their interaction probability in the model will be high. Hamoen et al. reported that ComK may bind to another ComK to form a homodimer. This is represented in the model using an interaction probability ρ (Table 2), with the following rule executed whenever a ComK finds itself next to another ComK agent, listed as neighbor here:
if neighbor = ComK then
random = generate random number between 0 and 1.
if random < ρ then
neighbor now moves with ComK
Initial agent quantities
Starting number at t0
If the probability threshold is not met, then the rule is not executed. The process is repeated at each time step, until eventually the probability threshold is met and the rule is executed.
The parameter space is quite large for this model and is divided into three types of parameter estimates: grid environment size, initial number of agents, and rule interaction probabilities. The parameters were estimated using a random sweep parameter estimation technique, where: (1) initial parameter estimates were made based on biological insight (see next paragraph), (2) repeated simulations were then run where parameter values were randomly changed, and (3) final parameter values were determined from simulations which best mimicked experimentally obtained results, as reported by Maamar et al..
Rule probabilities were initially estimated based on known interactions and then fit as simulations were run. For instance, if two molecules had a high affinity, then a high probability of interaction (e.g. 0.8) was estimated. For two molecules with a low affinity, a low probability (e.g. 0.2) was estimated. As simulations were run, rule probabilities were increased or decreased to speed up or slow down, respectively, a particular molecular reaction.
Intracellular agent-based model
Interacting agents, rule execution probabilities, and behavior
Repressor + Promoter
DegU + Promoter
ComX + Promoter
ComK + Promoter
ComK + Promoter + DegU
ComK + ComK
Ribosome + mRNA
ClpC/ClpP + MecA-ComK
ClpC/ClpP + MecA + ComK
ClpC/ClpP + MecA-ComS
ClpC/ClpP + MecA + ComS
ClpC/ClpP + MecA
ClpC/ClpP + MecA + ComK/ComS
Remove ComK or
New Cell if not starving
Remove Repressor if starving
ComK, ComS, ComX, DegU, MecA, ClpC/ClpP, mRNA, Ribosome, Repressor
Random walk, see text
Promoter + ComK dimer
Promoter + ComK tertramer
Promoter + ComX
Ribosome + mRNA
New ComK or
ComK transcription and translation
ComK binds at its own promoter as a tetramer, acting as its own transcription factor[6, 26]. Random transcription of comK is a key factor in the build-up of large amounts of ComK, which triggers transcription of the DNA transport (competence) genes[2, 27]. DegU binds to the comK promoter and strongly stimulates binding of ComK dimers to the comK promoter[28, 29]. More specifically, DegU binds in between the two ComK dimer binding sites and may possibly facilitate tetramerization of ComK at the comK promoter site by partial unwinding and bending of the DNA helix. Transcription can also occur in the absence of ComK. It is estimated that a cell exhibiting competence has on average 50,000 ComK dimers during stationary phase.
To simulate comK transcription in the intracellular model, agents representing ComK, the comK promoter, and DegU are modeled. Binding of ComK agents to the comK promoter agent is contingent upon what agents, if any, are already bound; in the model, this is reflected in the creation of binding rule probabilities (Table 3) calculated during the parameter estimation process. If no agents are bound to the promoter agent at the time a ComK dimer encounters it, the lowest binding probability is used. If DegU is already bound to the promoter agent, a higher probability of binding is used. Finally, the highest probability is used for binding the second ComK dimer if DegU and another ComK dimer are present.
Transcription probabilities follow in a similar fashion (Table 3). Transcription occurs at a very low probability when no ComK is bound – its probability increases with the addition of one or two bound ComK molecules. Activators and repressors (described in comK Transcriptional Regulation) disassociate upon successful completion of transcription. After the transcription rule is executed, an mRNA agent is generated, which will persist until it randomly degrades as defined by its “death” rule. Since the comK transcript has a strong Shine-Delgarno ribosome initiation sequence, the binding and translation probability of the ribosome agent is very high if it encounters a comK transcript (Table 3). Therefore, the presence of a single transcript is often enough to lead to the production of several ComK proteins, as reported in Maamar et al..
ComK transcriptional regulation
At present, there are three known comK transcriptional repressors: Rok, AbrB and CodY. The comK promoter site allows simultaneous binding of AbrB and ComK. The presence of AbrB acts to prevent binding of RNA polymerase, as does CodY. In the model, all three repressors are represented by a generic repressor agent that binds to the comK promoter agent. Nutrient limiting conditions down-regulate both AbrB and CodY, thus repressor agents are removed from the model when the cell agent reaches a starvation condition described below. Successful transcription and translation events eventually lead to the build-up of ComK, triggering the feedback loop shown in Figure 1[2, 27].
ComS transcription and regulation
ComS is a protein produced in response to quorum sensing (cell density)[33, 34]. Transcription of comS occurs in response to the quorum-sensing signaling pathway initiated by the ComX peptide. ComX is produced by the cell at a constant rate during growth and accumulates in the cell medium, reflecting cell density[1, 35]. For the purposes of the cell model, ComX is the post-translationally modified and cleaved extracellular end product absorbed by the cell. ComX initiates the activation of several proteins, which in turn initiate transcription of comS. The model represents this process by utilizing only a ComX agent. When the ComX agent binds to the comS promoter, comS is activated. The interaction probability given in Table 3 controls this rule, which represents a simplified version of the actual activation pathway. In reality, ComX is not the actual comS transcription factor. In addition, regulation of the quorum-sensing pathway is modeled by assuming a repressor agent acts at the comS promoter site, as shown in Table 2.
ComK and ComS post-translational regulation
The MecA/ClpC/ClpP protease complex degrades both ComK and ComS proteins. There are approximately 300 MecA/ClpC/ClpP molecules in a cell. MecA, an adapter protein, binds with either ComK or ComS, targeting the proteins for degradation by ClpC/ClpP. ComS competes with ComK for binding with MecA, with ComS having a higher affinity than ComK. If ComK is bound to MecA upon encountering ComS, ComK disassociates, targeting ComS for degradation instead. Because ComK is positively auto-regulated, protection from degradation by ComS results in an explosive increase in ComK synthesis[36, 37]. In this fashion, the up-regulation of ComS due to quorum sensing leads to an increased accumulation of ComK, and the cell transitions to the competence state. Thus, in the model, when ClpP/ClpC + MecA + ComK-bound agents encounters a ComS agent, ComK is released and ComS is bound.
This system was represented by implementing the interaction probabilities shown in Table 3. Since ComS has a stronger affinity to the adapter protein MecA, a higher binding probability is used for association with ComK during execution of the binding rule.
Intracellular model rules
The move rule simulates a random walk of an agent through the 3-D landscape to imitate molecular diffusion and Brownian motion. Each successful execution of the rule results in a one-step move of an agent to a randomly selected neighboring cell on the grid. In the 3-D grid of the intracellular model, there are 26 possible adjacent neighboring cells into which an agent could attempt to move. If the randomly selected cell is occupied, then the agent remains in place.
This rule is used for molecules that bind with other molecules when encountering another in a neighboring cell during the course of the random walk. For each agent, the bind rule searches for other agents at adjacent grid positions. When there is an adjacent agent for which a binding rule is defined, the two bind by moving together if the rule probability is met. An exception occurs when a MecA agent already bound to ComK encounters ComS – in this instance, ComK dissociates in favor of ComS. Agents that bind with one another are provided in Table 3 along with their associated rule probabilities determined during the parameter estimation process.
The transcription rule is executed by the promoter agent and results in the production of mRNA agents. The success of this rule depends on what is already bound to the promoter agent, to match the biology described above and as specified in Table 3.
The ribosome agent executes the translation rule to generate ComK agents or ComS agents depending on the type of the mRNA agent to which it is bound (Table 3).
Both the mRNA and ClpP/ClpC agents implement a death rule. For the mRNA agent, the death rule represents the random degradation of mRNA that occurs in the cell and is executed when unbound. In the case of the ClpP/ClpC protease complex, the death rule initiates the removal/death of the bound ComK or ComS agent. Table 3 lists the rule probabilities.
Culture agent-based model
In the cell culture model, agents interact with one another and their environment in a 2-D grid of 40 × 40 cells. Rules are executed in a random order for each iteration of the model, and probabilities calculated during the parameter estimation process determine whether or not a rule executed, as shown in Table 3.
Cell culture agents
There is technically only one agent type in the cell culture model, a cell ABM whose internal workings are described above. In the culture model, the cell ABM acts as an agent with the rules specified in Table 1. It interacts with its environment by consuming nutrients, and with other cell agents through production and consumption of the ComX pheromone. Nutrients and the ComX peptide are modeled by diffusion equations due to their high concentrations. Due to the way Repast Simphony is structured, the culture plate itself acts as a layer of immobile agents to allow for the execution of the diffusion rule. Consumption of either nutrients or ComX peptides by the cell agent are fed to the internal cell models, leading to the reduction of repressor agents or the increase of ComX agents, respectively. When more than 20 ComK agents were generated, the cell was considered to have reached competence (shown as a change in color in Figure 2). This threshold is an estimated parameter and was set to this small value to limit the consumption of computational resources.
Cell growth equation
The cell agent implements a growth function to control cell growth, division, and death. The growth function is based on the Logistic Map function: m n+1 = μm n -μm 2 n /k, where μ = 0.0058 is the growth rate, m represents the energy of the cell, and k is the maximum energy. m n is the value of the energy function at iteration n. In the culture ABM, m 0 was initialized to 5 and k was set to 16. The μm n term of the equation signifies an energy gain that occurs when consuming nutrients. The μm 2 n /k term signifies a decrease in energy, as it is assumed that basic metabolism within the cell consumes energy. The calculation of the growth function is split across two rules that are described further below: move and consumeNutrients. The life rule uses the value of the equation to determine whether the cell had accumulated enough energy to divide or not. In the move rule, energy is reduced by μm 2 n /k. If there are no nutrients at the agent’s location, there is no energy increase.
Once nutrients reach a low level (<1) such that cell agents could no longer consume them, the energy level steadily decreases instead of increasing via the move rule. At this stage, the cell growth equation is altered and energy is reduced by d/(k/2), where d is the death rate, d = 0.002.
Cell culture rules
The cell culture executes essentially one global rule, the diffusion rule. This rule executes the Repast Simphony diffusion algorithm on both the nutrient and peptide value layers, based on Rucker’s diffusion equation for Cellular Automata. For each cell in the grid, the difference between the current cell value and the weighted average of neighboring cells is calculated. The resulting value is then multiplied by the diffusion constant and added to the current cell value, thus ensuring that the concentrations within the grid cells move down the concentration gradient. The diffusion constant determined during the parameter estimation process is 0.1.
The cell agent executes a move rule that simulates bacterial chemotaxis so the cell agent moves towards the most favorable nutrient conditions. If nutrients are available, the cell agent remains at its location. If not, the cell agent randomly selects a free neighboring location with a higher concentration of nutrients. The energy reduction portion of the cell growth equation is implemented when this rule is executed due to the assumption that there is a metabolic cost to a cell’s energy at each iteration.
The move rule also determines when nutrients can no longer be consumed at the cell agent’s current location, and thus disables repressor agents within the cell ABM. A randomly selected repressor agent is removed from the model if a probability threshold (0.001) determined during the parameter estimation process is met.
The consumePeptide and consumeNutrients rules cause the consumption of one molecule from the grid location of the cell agent. When a ComX peptide is consumed, it is added as an agent to the cell ABM. When a nutrient is consumed, a gain in energy occurs, as described above. The consumeNutrient rule is executed at every iteration, while the consumePeptide rule is executed every 50 iterations, if a probability of 0.8 is met.
The generatePeptide rule produces one molecule of the ComX peptide, which is added to the concentration at the current grid location. This rule is an approximation of the constant production of the ComX peptide described in the literature. Unlike the other rules, the generatePeptide rule is executed every 100 iterations and a ComX peptide is generated with a probability of 0.8. This number of iterations was selected during the parameter estimation process to ensure a gradual production of the peptide. The ComX peptides produced in this manner are independent of the ComX agents residing within the cell ABM agent.
The life rule determines whether the cell is ready to divide when the energy exceeds a predefined threshold of 15, which is one fewer than k, the maximum energy (see growth equation parameters above). Cells that exhibit the competence state do not divide. During the division process, a new cell ABM is created and added to the grid in a neighboring, adjacent grid location. The daughter cell receives half the energy of the parent cell, thereby reducing the parent’s energy by half. The daughter cell will follow intracellular model startup and agent initialization as described in previous sections and in Table 2. However, some of the initial agent quantities are handled differently.
To model inheritance, an arbitrary plane is randomly chosen which bisects the parent cell ABM through its center. ComK, ComS, ComX and mRNA agents ‘above’ the plane remain in the parent cell and agents ‘below’ the plane become a part of the daughter cell. The daughter cell is then placed in a randomly selected grid location adjacent to the parent cell. If that location is occupied, then the shove rule is executed. The simulation assures that the associated genes (promoter agents) are equally distributed between mother and daughter to represent the non-stochastic genome segregation in real cells.
The life rule is also responsible for determining the death of a cell agent. The effect of this rule is shown during the death phase of the bacterial growth curve. If the cell agent’s energy becomes very low (< 0.5) and the probability threshold of 0.0001 is met, then the cell agent ‘dies’ and is removed from the model.
The shove rule is intended to displace cell agents one step to an adjacent, randomly selected, neighboring position if more than one cell agent occupies its current location. As each agent executes this rule, a one step displacement of cell agents ripples through a group of adjacent cell agents until there is room for all cell agents in the culture model grid.
At model start-up, agents were placed at random locations within the simulated 3-D grid environment, Figure 3c. Several simulations were run to assess outcomes given different initial configurations.
Cell culture level modeling
One goal of the study was to leverage the intracellular model into a multi-scale model of cell culture level interactions, to examine the effects of signaling and competition for nutrients on the outcome of competence. The cell-level model consisted of whole cell agents representing a cell’s interaction with external environmental factors such as nutrients and the quorum-sensing pheromone ComX. Agents were placed randomly in a 2-D grid environment (Figure 3b), much like the random placement of agents in the 3-D intracellular tier (Figure 3c).
In this multi-scale model, external environmental conditions influenced the intracellular conditions and cell-level outcomes, and intracellular conditions fed back upon the environment and other agents within it. For all multi-scale simulations, the initial intracellular model quantities of ComK, ComS, and ComX agents were randomly determined from thresholds as explained in Model Implementation. There were no comK mRNA agents placed at the start of a model run, but these agents were randomly generated via transcription during a simulation.
Successful execution of transcription and translation rules eventually led to the build-up of ComK, triggering the feedback loop shown in Figure 1[2, 27]. As competition between ComK and ComS proteins is vital in determining the competence state, their production was monitored. Considerably fewer molecules were used in the simulations compared to actual biology due to limitations in computational resources. As a result, a cell was considered to be in the competence state when it had more than 20 ComK proteins – a value determined during the parameter search stage that met the 10–20% competent cells criterion. In reality, approximately 50,000 ComK dimers have been observed in a competent cell.
Competence outcome at the intracellular level: the impact of random spatial-temporal agent arrangement
Initial simulations of the intracellular model were run repeatedly outside of the full model environment to test if spatial arrangement of molecules might contribute to the resulting competence state. Intracellular model simulations were run with identical parameters and concentrations of agents (Table 2, maximum values for ComK, ComS, ComX and DegU were used in these simulations), the sole difference between each run being the random spatial placement of the agents before the simulation started. We ran the simulation 100 times to obtain statistics regarding the overall rate of competence emergence.
We repeatedly observed that simulations with identical initial parameters and population sizes resulted in significantly different outcomes, suggesting that the initial arrangement of molecules indeed has a strong impact on competence expression. Given the parameter settings, the over-production of ComK was observed in approximately 3–5 out of 100 simulations (Figure 5b) and the production of ComK outpaced repressor protein activity. Though 3–5% is a lower rate of competence expression than the 10–20% that is typically seen in the bacteria, the model was run without explicit simulation of cell growth, division, or starvation conditions that influence a higher density of competence emergence.
Another outcome of the simulations was that the chance encounter of a single comK mRNA transcript with a ribosome before it was degraded would produce a ComK protein. This agrees with prior in vivo results indicating that an increased probability of competence switching can be driven by fluctuation of the ComK protein by a few ComK mRNA molecules. On average, one mRNA was observed at any given time in the model comparied to > 20 during growth phase and on average 1 mRNA during stationary phase in vivo. So cell fate in the ABM was determined simply by the location of the molecule and its chance encounter with a ribosome. These results provide evidence that a random, spatial arrangement of molecules may be a major contributor to the variable comK gene expression (noise) underlying the competence phenotype, and could lead to experiments that further decode the enigma of competence expression.
Competence outcome at the cell culture level: the impact of nutrient limitation
It has been previously shown that nutrient limitation and cell culture density increase the propensity of B. subtilis cells to enter the competence state. In our intracellular model, only a small fraction of the simulations ended in the competence state since the model did not include limitations comparable to real-life conditions. However, within actual cell cultures, once stationary phase is reached and nutrients become limited, a much higher fraction of cells emerge with the competence phenotype - as many as 20%. We sought to determine what bottom-up assumptions drove this number by running simulations of the complete model – intracellular ABMs acting as agents in the culture model.
At each simulation run, the culture model was seeded with 20 intracellular ABMs placed in random locations in the 2-D grid environment. Each of the 20 initial intracellular models began with quantities of ComK, ComS, ComX, and DegU determined by a random draw from a uniform distribution from within a defined threshold, and with 0 mRNA agents. Cells grew and divided, with division resulting in a random partitioning of molecular contents among the progeny. All cells were tracked throughout the simulation, and progeny progressed from 20 initial cells to a maximum of 877 cells in one simulation, with 137 (15.6%) showing competence at simulation termination.
Since individual simulations would often result in distinct outcomes, we ran the model repeatedly to obtain average statistics for competence-switching behavior. In six simulations, the initial seed of 20 randomly placed cells produced an average of 867 cells. We limited the available “plate size” and nutrient concentration for culture growth to limit the computing to feasible time spans. For each of the cells produced, a complete intracellular model was running, which meant that a full simulation running on a fast desktop computer took approximately two months or more to complete. Improved parallelism will reduce run times in the future.
Additional file 1: Movie of Cell Culture Model. This movie shows the cell culture ABM initially seeded with 20 cells (blue). Green indicates nutrient-rich cells. Cells grow and divide as they consume nutrients. The formerly nutrient-rich cells darken as nutrient concentrations decrease. Cells move following the nutrient gradient as nutrients are depleted. When cells die, nutrients are released into the environment. The movie halts during stationary phase. (MOV 1 MB)
Competence outcome at the cell culture level: epigenetic heritability
Nutrient limitation and cell density does not, however, explain why those particular 10–20% of cells that express competence do so. Following from the result showing that initial spatial distribution of molecules may affect the outcome of the intracellular model, we postulated that some level of epigenetic heritability might exist in competence switching. The rationale is that if a parent cell has an increased quantity of ComK compared to the average, it might be expected that as the cell divides, the resulting progeny may also have elevated ComK levels. Evidence for this type of heritability in the biological system was presented by Veening et al.. The precursor to competence – accumulation of ComK – may lead to an earlier switch than would occur by random chance alone, a result which would be clearly identifiable in our cell culture model.
Our interest in phenotype switching derives from studies of antimicrobial-tolerant bacterial phenotypes. These are non-inherited phenotypes that confer tolerance towards many known antibiotic drugs. In bacterial populations, cells exhibiting tolerant phenotypes exist as a small fraction of the population, with the quantity determined in part by growth conditions. For example, stationary phase growth induces an increase in the fraction of phenotypically drug-tolerant cells. However, at the present time, little is known about the mechanisms underlying antimicrobial tolerance in bacteria. So, to explore mechanisms that could explain variable gene expression and to further our knowledge of phenotypic bistable switching, we turned to competence switching in B. subtilis, where the key molecular players are known and well-studied.
Our B. subtilis ABM is a reduced scale model tuned to reflect the bacteria biology and the competence phenotype. Through the Cell Growth Equation (see Model Implementation section) the model closely approximates bacterial growth curves. In addition, cells only turned competent in bulk once stationary phase was entered. Although quantities of molecular agents modeled were a small fraction of the molecules in a live system, the model was able to approximate a real system. We have found similar results in a separate effort modeling chemotaxis in E. coli with an ABM, where preserving biological ratios of molecules was more important than preserving absolute quantities to produce biologically realistic results. However, in the case of this model we found that chance, random interactions between agents had the most effect on the behavior of the model.
In our model, spatial-temporal interactions – not molecular quantities or agent rule execution probabilities – may be the key contributor to the emergence of the competence phenotype, as transcription and translation occur due to chance encounters by the molecules involved. Though it is not yet possible to verify these simulations in vivo, the biological system likely shares similar properties where the apparent randomness of the competence transitions is derived directly from the randomness of spatial distribution in competence-related molecules. Noise in a biological system is typically thought of as a series of events, but our results show it is possible that random spatial arrangement of competence-determining molecules may be its major contributing factor.
In addition to random spatial arrangements, random encounters between key agents also appeared to be important in the emergence of competence. It is interesting to note that changing the quantities of the initial number of agents in the model had more of an effect on the model outcome than the rule probabilities themselves. If there are greater numbers of an agent, then there is a higher chance for encounters between agents to occur. During the parameter estimation process, it was observed that changes in the probability values did not greatly impact model outcomes, although higher probabilities did speed up certain reactions (e.g. translation), which resulted in the production of more agents. This seems logical considering that the rule probabilities do not come into effect until two agents encounter one another, which can occur more frequently as the quantities of agents increase. If there are a large number of agents, there are more chances for them to encounter one another. As parameters have been reported to lack sensitivity in many biological models, our parameterization method is somewhat un-formalized as we focused more on looking for biological interpretations of our model after it had reached a life-like representation of the competence switch.
The cell division trees that monitored quantities of ComK displayed an obvious pattern of inheritance. Cells that inherited higher amounts of ComK transcripts and proteins from the parent cell tended to, in turn, pass on higher quantities to their children cells. However, the model also suggests what may be a basic but important facet of the emergence of competence: that dilution of competence-determining molecules during cell division may act to regulate the emergence of competence. Except for the lineages that switched to the competent state early due to inheriting elevated levels of the ComK protein, the remaining competent cells only became so after nutrient limitation and repressor agents inhibited cell division long enough for sufficient quantities of ComK agents to accumulate. The quantities of molecules inherited in this pathway acted in an epigenetic fashion upon subsequent phenotypic outcomes. This explains why competence in B. subtilis typically emerges during stationary phase when nutrients are limited, due to build-up of elements like ComK in the competence gene regulatory network.
The model may therefore reveal an important feature of how the growth stage regulates competence by diluting the mRNA and proteins essential to switching to the competent state. Although we have not verified this experimentally, Roostalu et al. reported fluctuations in a reporter protein after cell division that supports the theory that cell division acts as a dilution event for the competence phenotype. During exponential growth phase, not enough time would pass between cell divisions to allow ComK to build up to sufficient levels, so competence would not occur. While we cannot know whether this is the complete explanation for how competence is limited to stationary phase in B. subtilis cultures, it appears to be a sufficient explanation to guide further experiments.
This study suffers from the same limitations as any model-based research, in that all models are only representations of a biological system. Nevertheless, the design of the model facilitates a thought experiment to represent and clarify the workings behind very complex cell behavior. This discrete model of the competence phenotype provides a readily comprehensible view into each cell’s behavior and provides the ability to monitor the variation of molecular concentrations involved in regulating competence. The resulting model is biologically intuitive, with ready translation from biological facts or hypotheses into the model and back. As we used this model to educate ourselves on the mechanisms of the competence phenotype, we see great potential in ABMs as educational tools in the future because they are straightforward to build, visualize, and comprehend.
We thank Wendy Spitzer and Ashley Lundqvist who provided technical writing services. We also wish to thank Jameson Miller and Miles Parker for sharing their agent-based modeling skill sets, as well as Dr. Kevin Ramkissoon and Dr. Hsun-Cheng Su for fruitful discussions about bacterial growth curves. We thank the members of the Repast-Interest mailing list for their support of the Repast Simphony platform.
A portion of this work was supported by F31DE019594 from the National Institute of Dental & Craniofacial Research to SMS. This work was also supported by NIH R01RR020823 to MCG. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institute of Dental & Craniofacial Research, the National Center for Research Resources, or the National Institutes of Health.
- Draskovic I, Dubnau D: Competence for genetic transformation. The dynamic bacterial genome. Volume 8. Edited by: Mullany P. 2005, Cambridge, UK; New York: Cambridge University Press, 235-273.View ArticleGoogle Scholar
- Maamar H, Raj A, Dubnau D: Noise in Gene Expression Determines Cell Fate in Bacillus subtilis. Science. 2007, 317 (5837): 526-529. 10.1126/science.1140818.View ArticlePubMedGoogle Scholar
- Henderson IR, Owen P, Nataro JP: Molecular switches: the ON and OFF of bacterial phase variation. Mol Microbiol. 1999, 33 (5): 919-932. 10.1046/j.1365-2958.1999.01555.x.View ArticlePubMedGoogle Scholar
- Dubnau D, Losick R: Bistability in bacteria. Mol Microbiol. 2006, 61 (3): 564-572. 10.1111/j.1365-2958.2006.05249.x.View ArticlePubMedGoogle Scholar
- Webb JS, Lau M, Kjelleberg S: Bacteriophage and Phenotypic Variation in Pseudomonas aeruginosa Biofilm Development. J Bacteriol. 2004, 186 (23): 8066-8073. 10.1128/JB.186.23.8066-8073.2004.PubMed CentralView ArticlePubMedGoogle Scholar
- van Sinderen D, Luttinger A, Kong L, Dubnau D, Venema G: Hamoen L: comK encodes the competence transcription factor, the key regulatory protein for competence development in Bacillus subtilis. Mol Microbiol. 1995, 15 (3): 455-462. 10.1111/j.1365-2958.1995.tb02259.x.View ArticlePubMedGoogle Scholar
- Hamoen LW, Venema G, Kuipers OP: Controlling competence in Bacillus subtilis: shared use of regulators. Microbiology. 2003, 149 (1): 9-17. 10.1099/mic.0.26003-0.View ArticlePubMedGoogle Scholar
- Maughan H, Nicholson WL: Stochastic Processes Influence Stationary-Phase Decisions in Bacillus subtilis. J Bacteriol. 2004, 186 (7): 2212-2214. 10.1128/JB.186.7.2212-2214.2004.PubMed CentralView ArticlePubMedGoogle Scholar
- Veening J, Smits W, Kuipers O: Bistability, Epigenetics, and Bet-Hedging in Bacteria. Annu Rev Microbiol. 2008, 62: 193-210. 10.1146/annurev.micro.62.081307.163002.View ArticlePubMedGoogle Scholar
- Drenkard E, Ausubel FM: Pseudomonas biofilm formation and antibiotic resistance are linked to phenotypic variation. Nature. 2002, 416 (6882): 740-743. 10.1038/416740a.View ArticlePubMedGoogle Scholar
- Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: from theories to phenotypes. Nat Rev Genet. 2005, 6 (6): 451-464. 10.1038/nrg1615.View ArticlePubMedGoogle Scholar
- Smits W, Kuipers O, Veening J: Phenotypic variation in bacteria: the role of feedback regulation. Nat Rev Micro. 2006, 4 (4): 259-271. 10.1038/nrmicro1381.View ArticleGoogle Scholar
- Süel G, Garcia-Ojalvo J, Liberman L, Elowitz M: An excitable gene regulatory circuit induces transient cellular differentiation. Nature. 2006, 440 (7083): 545-550. 10.1038/nature04588.View ArticlePubMedGoogle Scholar
- Süel G, Kulkarni RP, Dworkin J, Garcia-Ojalvo J, Elowitz M: Tunability and noise dependence in differentiation dynamics. Science. 2007, 315 (5819): 1716-1719. 10.1126/science.1137455.View ArticlePubMedGoogle Scholar
- Schultz D, Ben Jacob E, Onuchic J, Wolynes P: Molecular level stochastic model for competence cycles in Bacillus subtilis. Proc Natl Acad Sci USA. 2007, 104 (45): 17582-17587. 10.1073/pnas.0707965104.PubMed CentralView ArticlePubMedGoogle Scholar
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81 (25): 2340-2361. 10.1021/j100540a008.View ArticleGoogle Scholar
- Kreft JU, Picioreanu C, Wimpenny JW, van Loosdrecht MC: Individual-based modelling of biofilms. Microbiology. 2001, 147: 2897-2912.View ArticlePubMedGoogle Scholar
- Goryachev A, Toh D, Wee K, Lee T, Zhang H, Zhang L: Transition to quorum sensing in an Agrobacterium population: A stochastic model. PLoS Comput Biol. 2005, 1 (4): e37-10.1371/journal.pcbi.0010037.PubMed CentralView ArticlePubMedGoogle Scholar
- Xavier J, Foster K: Cooperation and conflict in microbial biofilms. Proc Natl Acad Sci USA. 2007, 104 (3): 876-881. 10.1073/pnas.0607651104.PubMed CentralView ArticlePubMedGoogle Scholar
- Johnson LR: Microcolony and biofilm formation as a survival strategy for bacteria. J Theor Biol. 2008, 251 (1): 24-34. 10.1016/j.jtbi.2007.10.039.View ArticlePubMedGoogle Scholar
- Nadell C, Xavier J, Levin S, Foster K: The evolution of quorum sensing in bacterial biofilms. Plos Biol. 2008, 6 (1): e14-10.1371/journal.pbio.0060014.PubMed CentralView ArticlePubMedGoogle Scholar
- D'Agata EM, Magal P, Olivier D, Ruan S, Webb GF: Modeling antibiotic resistance in hospitals: the impact of minimizing treatment duration. J Theor Biol. 2007, 249 (3): 487-499. 10.1016/j.jtbi.2007.08.011.PubMed CentralView ArticlePubMedGoogle Scholar
- Murphy JT, Walshe R, Devocelle M: A computational model of antibiotic-resistance mechanisms in methicillin-resistant Staphylococcus aureus (MRSA). J Theor Biol. 2008, 254 (2): 284-293. 10.1016/j.jtbi.2008.05.037.View ArticlePubMedGoogle Scholar
- Miller J, Parker M, Bourret RB, Giddings MC: An Agent-Based Model of Signal Transduction in Bacterial Chemotaxis. PLoS One. 2010, 5 (5): e9454-10.1371/journal.pone.0009454.PubMed CentralView ArticlePubMedGoogle Scholar
- Repast Suite: [http://repast.sourceforge.net/] 
- Hamoen LW, Van Werkhoven AF, Bijlsma JJE, Dubnau D, Venema G: The competence transcription factor of Bacillus subtilis recognizes short A/T-rich sequences arranged in a unique, flexible pattern along the DNA helix. Genes Dev. 1998, 12 (10): 1539-1550. 10.1101/gad.12.10.1539.PubMed CentralView ArticlePubMedGoogle Scholar
- Leisner M, Stingl K, Radler JO, Maier B: Basal expression rate of comK sets a ‘switching-window’ into the K-state of Bacillus subtilis. Mol Microbiol. 2007, 63 (6): 1806-1816. 10.1111/j.1365-2958.2007.05628.x.View ArticlePubMedGoogle Scholar
- Ogura M, Tanaka T: Bacillus subtilis DegU acts as a positive regulator for comK expression. FEBS Lett. 1996, 397 (2–3): 173-176.View ArticlePubMedGoogle Scholar
- Hamoen L, Van Werkhoven A, Venema G, Dubnau D: The pleiotropic response regulator DegU functions as a priming protein in competence development in Bacillus subtilis. Proc Natl Acad Sci USA. 2000, 97 (16): 9246-9251. 10.1073/pnas.160010597.PubMed CentralView ArticlePubMedGoogle Scholar
- Smits W, Hoa T, Hamoen L, Kuipers O, Dubnau D: Antirepression as a second mechanism of transcriptional activation by a minor groove binding protein. Mol Microbiol. 2007, 64 (2): 368-381. 10.1111/j.1365-2958.2007.05662.x.View ArticlePubMedGoogle Scholar
- Hamoen LW, Kausche D, Marahiel MA, van Sinderen D, Venema G, Serror P: The Bacillus subtilis transition state regulator AbrB binds to the −35 promoter region of comK. FEMS Microbiol Lett. 2003, 218 (2): 299-304. 10.1111/j.1574-6968.2003.tb11532.x.View ArticlePubMedGoogle Scholar
- Ratnayake-Lecamwasam M, Serror P, Wong K, Sonenshein AL: Bacillus subtilis CodY represses early-stationary-phase genes by sensing GTP levels. Genes Dev. 2001, 15 (9): 1093-1103. 10.1101/gad.874201.PubMed CentralView ArticlePubMedGoogle Scholar
- D'Souza C, Nakano MM, Zuber P: Identification of comS, a gene of the srfA operon that regulates the establishment of genetic competence in Bacillus subtilis. Proc Natl Acad Sci U S A. 1994, 91 (20): 9397-9401. 10.1073/pnas.91.20.9397.PubMed CentralView ArticlePubMedGoogle Scholar
- Hamoen LW, Eshuis H, Jongbloed J, Venema G, Sinderen D: A small gene, designated comS, located within the coding region of the fourth amino acid-activation domain of srfA, is required for competence development in Bacillus subtilis. Mol Microbiol. 1995, 15 (1): 55-63. 10.1111/j.1365-2958.1995.tb02220.x.View ArticlePubMedGoogle Scholar
- Schneider KB, Palmer TM, Grossman AD: Characterization of comQ and comX, Two Genes Required for Production of ComX Pheromone in Bacillus subtilis. J Bacteriol. 2002, 184 (2): 410-419. 10.1128/JB.184.2.410-419.2002.PubMed CentralView ArticleGoogle Scholar
- Turgay K, Hahn J, Burghoorn J, Dubnau D: Competence in Bacillus subtilis is controlled by regulated proteolysis of a transcription factor. EMBO J. 1998, 17 (22): 6730-6738.PubMed CentralView ArticlePubMedGoogle Scholar
- Prepiak P, Dubnau D: A Peptide Signal for Adapter Protein-Mediated Degradation by the AAA + Protease ClpCP. Mol Cell. 2007, 26 (5): 639-647. 10.1016/j.molcel.2007.05.011.PubMed CentralView ArticlePubMedGoogle Scholar
- Burdett ID, Kirkwood TB, Whalley JB: Growth kinetics of individual Bacillus subtilis cells and correlation with nucleoid extension. J Bacteriol. 1986, 167 (1): 219-230.PubMed CentralPubMedGoogle Scholar
- Fall CP, Marland ES, Wagner JM, Tyson JJ: Computational Cell Biology. 2004, New York, NY: Springer-Verlag New York, Inc.View ArticleGoogle Scholar
- Continuous Valued Cellular Automata in Two Dimensions: [http://www.cs.sjsu.edu/faculty/rucker/capow/santafe.html] 
- Shapiro L, McAdams HH, Losick R: Why and how bacteria localize proteins. Science. 2009, 326 (5957): 1225-1228. 10.1126/science.1175685.View ArticlePubMedGoogle Scholar
- Maamar H, Dubnau D: Bistability in the Bacillus subtilis K-state (competence) system requires a positive feedback loop. Mol Microbiol. 2005, 56 (3): 615-624. 10.1111/j.1365-2958.2005.04592.x.View ArticlePubMedGoogle Scholar
- Hoa TT, Tortosa P, Albano M, Dubnau D: Rok (YkuW) regulates genetic competence in Bacillus subtilis by directly repressing comK. Mol Microbiol. 2002, 43 (1): 15-26. 10.1046/j.1365-2958.2002.02727.x.View ArticlePubMedGoogle Scholar
- Veening J, Stewart EJ, Berngruber TW, Taddei F, Kuipers OP, Hamoen LW: Bet-hedging and epigenetic inheritance in bacterial cell development. Proc Natl Acad Sci USA. 2008, 105 (11): 4393-4398. 10.1073/pnas.0700463105.PubMed CentralView ArticlePubMedGoogle Scholar
- R Development Core Team: R: A Language and Environment for Statistical Computing. [http://www.r-project.org] 
- Levin B, Rozen D: Non-inherited antibiotic resistance. Nat Rev Micro. 2006, 4 (7): 556-562. 10.1038/nrmicro1445.View ArticleGoogle Scholar
- Keren I, Kaldalu N, Spoering A, Wang Y, Lewis K: Persister cells and tolerance to antimicrobials. FEMS Microbiol Lett. 2004, 230 (1): 13-18. 10.1016/S0378-1097(03)00856-5.View ArticlePubMedGoogle Scholar
- Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally sloppy parameter sensitivities in systems biology models. PLoS Comput Biol. 2007, 3 (10): 1871-1878.PubMedGoogle Scholar
- Roostalu J, Joers A, Luidalepp H, Kaldalu N, Tenson T: Cell division in Escherichia coli cultures monitored at single cell resolution. BMC Microbiol. 2008, 8 (1): 68-10.1186/1471-2180-8-68.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.