article: Simulation of cellular behavior on multiple scales requires models that are sufficiently detailed to capture central intracellular processes but at the same time enable the simulation of entire cell populations in a computationally cheap way. In this paper we present RapidCell, a hybrid model of chemotactic Escherichia coli that combines the Monod-Wyman-Changeux signal processing by mixed chemoreceptor clusters, the adaptation dynamics described by ordinary differential equations, and a detailed model of cell tumbling. Our model dramatically reduces computational costs and allows the highly efficient simulation of E. coli chemotaxis. We use the model to investigate chemotaxis in different gradients, and suggest a new, constant-activity type of gradient to systematically study chemotactic behavior of virtual bacteria. Using the unique properties of this gradient, we show that optimal chemotaxis is observed in a narrow range of CheA kinase activity, where concentration of the response regulator CheY-P falls into the operating range of flagellar motors. Our simulations also confirm that the CheB phosphorylation feedback improves chemotactic efficiency by shifting the average CheY-P concentration to fit the motor operating range. Our results suggest that in liquid media the variability in adaptation times among cells may be evolutionary favorable to ensure coexistence of subpopulations that will be optimally tactic in different gradients. However, in a porous medium (agar) such variability appears to be less important, because agar structure poses mainly negative selection against subpopulations with low levels of adaptation enzymes. RapidCell is available from the authors upon request.
  One of the central questions of modern systems biology is the influence of microscopic parameters of a single cell on the behavior of a cell population, a common problem in multi-scale modeling. In terms of bacterial chemotaxis, this issue can be formulated as the influence of signaling network parameters on the spatiotemporal dynamics of a population in various gradients of chemoattractants. The problem of efficient multi-scale simulation imposes strict requirements on the model: it should be maximally detailed to grasp the main features of the signaling network yet computationally cheap to simulate large numbers of bacteria. Chemotaxis plays an important role in microbial population dynamics. Chemotactic bacteria in a nonmixed environment—that is in presence of nutrient gradients—have significant growth advantages, as shown experimentally for different bacterial species [1]–[4]. Modeling of microbial population dynamics indicates that motility and chemotactic ability can be as important for evolutionary competition as cell growth rate [5], [6]. Escherichia coli is an ideal organism for chemotaxis modeling, because of the rich experimental information collected over years of extensive research. In common with many other bacteria, E. coli can migrate towards high concentrations of attractants and away from repellents. In the adapted state, cells perform a random walk, which becomes biased in the presence of a spatial gradient. This swimming bias is based on temporal comparisons of attractant concentrations during cell runs. If the direction of a run is favorable, i. e. up the attractant gradient or down the repellent gradient, the run become longer. Between runs, the cell tumbles and reorients for the next run [7]. Chemotaxis in E. coli is mediated by an atypical two-component signal transduction pathway (for recent reviews see [8], [9]). Ligand molecules bind to clusters of transmembrane receptors, which are in complex with the histidine kinase CheA and the adaptor CheW. Each receptor can be either active or inactive, depending on ligand binding and the methylation level. The active receptor activates CheA, eliciting downstream phosphorylation of the response regulator CheY. Phosphorylated CheY (CheY-P) is dephosphorylated by CheZ. Receptors can be methylated by the methyltransferase CheR and demethylated by the methylesterase CheB, and methylation regulates the receptor activity. The methylation of receptors provides a sort of chemical ‘memory’, which allows the cell to compare the current ligand concentration with the past. Phosphorylation of CheB by CheA provides a negative feedback to the system, although it appears nonessential for exact adaptation [10], [11]. Phosphorylated molecules of CheY-P freely diffuse through the cytoplasm and bind to the flagellar motor protein FliM, causing motors to switch from CCW to CW rotation. Switching of the motors to the CW state results in a tumble and reorientation, whereas the CCW rotation corresponds to straight runs. A number of mathematical models of chemotaxis have been proposed [10], [12]–[18], including two recent programs that simulate cell motion along with the intracellular pathway dynamics: AgentCell [19], which is based on the StochSim pathway simulator [20]–[22], and E. solo [23], which is based on the BCT simulator [24]–[26]. The current version of AgentCell (2. 0) simulates the whole pathway stochastically, making it thus computationally very expensive. The E. solo program simulates the pathway using about 90 ordinary differential equations (ODEs). However, simulation of large bacterial populations on long time scales requires computationally cheaper models. It was recently shown using fluorescence resonance energy transfer (FRET) that the amplitude of the initial CheY-P response can be described by a Hill function of a relative change in receptor occupancy during stepwise ligand stimulation [27]. Recent modeling efforts [12], [28], [29] show that a mixed-cluster Monod-Wyman-Changeux (MWC) model of strongly coupled receptors is consistent with the FRET data, and can account for the observed sensitivity and precise adaptation over a wide range of ligand concentrations. The amplitude of pathway excitation can therefore be determined using several algebraic equations describing the free energy of the cluster. In our model (Figure 1A), we employed the MWC model for a mixed receptor cluster [12] with a mean-field approximation for adaptation kinetics [30]. Due to its hybrid approach, the model allowed us to reduce the computational costs dramatically, while keeping the main quantitative characteristics of the cell response (methylation level, relative CheY-P concentration, motor bias) consistent with experimental data. To couple the bias of individual motors to the probability of tumbling, we applied a voting model for several independent motors, based on detailed experimental investigation of tumbling mechanics [31]. These components were combined into a new simulator for E. coli chemotaxis—RapidCell, which uses a hybrid pathway simulation instead of a fully stochastic or ODE approach, and is therefore computationally cheap. This allows the simulation of populations of 104–105 cells on a time scale of hours using a desktop computer. To study the dependence of chemotaxis on gradient strength in a systematic way, we propose a new—constant-activity—gradient which ensures a constant average CheY-P level and cellular drift velocity along the gradient, in contrast to commonly used Gaussian and linear gradients. We show that the MWC model gives an approximately constant response over a wide range of ligand concentrations. Though purely theoretical, such a gradient serves as a perfect in silico assay to study the chemotactic properties of cells. The chemotaxis pathway is robust to changes in network parameters and intracellular protein concentrations [10], [15], [32]. This enables efficient chemotaxis with varying levels of intracellular components and under perturbations from extracellular environment. However, adaptation time is not robust [10], [11], [33], [34] and varies even among genetically identical cells in a population because of stochastic variations in gene expression and low copy numbers of the adaptation enzymes. Our simulations predict that in liquid media for any given gradient steepness, there is an optimal adaptation rate that provides the highest cellular drift velocity. We suggest a simple mechanism for this phenomenon: the optimal rate of adaptation is observed in a narrow range of kinase activity, where the average CheY-P level fits the operating range of the flagellar motor. In this range, the relation between CheY-P and motor bias is approximately linear, and cells perform chemotaxis with the highest efficiency. The situation is different for cells swimming in agar. Here, the optimal range of motor bias appears to be very narrow and just slightly higher than in the non-stimulated state. Due to the porous structure of agar, cells with a higher CCW motor bias stay trapped for a longer time, thus negating advantage in chemotactic efficiency. This leads to a strong selection against cells which adapt slowly and therefore tend to overreact to chemotactic stimulation. On the other hand, chemotaxis in agar poses only a weak selection against cells with a high adaptation rate. Our simulations suggest that in liquid media the variability in protein levels among cells may be advantageous for bacterial populations on a long time scales. In a nonmixed environment with different food sources and gradient intensities, such variability can help the whole population to respond to different gradients more readily, due to positive selection of subpopulations with optimal levels of adaptation enzymes in a given gradient.
  During a run, the cell is assumed to move with a constant speed v = 20 µm/s, while the direction of motion is affected by rotational diffusion [7], [42]. After each time step, the running direction is changed by adding a stochastic component with normal distribution and diffusion coefficient Dr = 0. 062 rad2 s−1 [42]. The virtual cells are swimming in a 2D environment with a predefined attractant concentration field S (x, y, t). The domain geometry is either rectangular or circular, with reflecting walls. The simulation time was chosen to be short enough to avoid boundary effects. The rectangular domain is within (0, xmax) × (0, ymax), and the circular domain within (0, rmax), with xmax = ymax = 2rmax = 20 mm. The output file of the RapidCell program contains the key characteristics of the intracellular state (CheY-P level, methylation state, motor bias) and the geometric characteristics of cell motion (position and orientation). The model was implemented using Java classes similar to AgentCell [19], but with simplified architecture. The algorithm is implemented as a discrete-time Monte Carlo scheme with time step Δt = 0. 01 s. For random-number generation, we used external Java libraries [49], [50]. The code was written using Eclipse SDK (www. eclipse. org). The output data were analyzed with MATLAB (The MathWorks, MA).
  To test our model (Figure 1A), we compared cellular behavior in the proposed universal constant-activity gradient with other gradients, observing the single cell swimming (Figure 1B) and the behavior of large populations. The key characteristics we consider are the CheY-P concentration and the drift velocity along the gradient. We used the constant-activity gradient to study the effect of adaptation rate on chemotactic efficiency. For this purpose, we simulated homogeneous populations consisting of cells with the same adaptation rate. In a fixed constant-activity gradient, the population drift velocity depends on adaptation rate in a unimodal manner (Figure 5A). A zero level of adaptation enzymes (non-adapting cells) results in a low drift velocity, though it is clearly distinguishable from non-chemotactic behavior. A proportional increase of adaptation rate improves cellular drift velocity up to a certain maximum, after which it slowly declines again. Extremely high adaptation rates, more than 100 times higher than wild-type, make the cells non-chemotactic (Figure 5A). To study chemotactic efficiency as a function of gradient steepness, cells were simulated in six constant-activity gradients with the steepness changing 64-fold, from 1. 14 to 72. 88×10−3, (Figure 5B). In each gradient, we determined the optimal adaptation rate, at which cellular drift velocity reaches its maximum. The simulated drift velocities are in the same range as those measured experimentally for E. coli in steep gradients (7 µm s−1) [53]. Our simulations indicate that experimental cell-drift velocities are inlikely to exceed 15 µm s−1, corresponding to an extremely steep and short-scale gradient. In very weak gradients, the drift velocity can be as low as 2. 5 µm s−1, still distinguishable from the non-chemotactic cellular drift (0. 8 µm s−1). Interestingly, we observed that the optimal adaptation rate rises in proportion with the gradient steepness (Figure 5B). To investigate the latter effect in more detail, we varied the adaptation rate from 0 to 10-fold relative to the wild-type. In steeper gradients, the optimal adaptation rate is indeed higher (Figure 6A), and the peak of the drift velocity becomes less sharp. To find the reason for the observed dependence between the gradient steepness and optimal adaptation rate, we tracked the average CheY phosphorylation levels of the virtual cells. As one can see in Figure 6A and 6B, in all gradients the 90%-intervals around the velocity peaks correspond to adaptation rate intervals [0. 1,0. 5], [0. 4,1. 5], [1,3], respectively. These three intervals fall into to the same interval [0. 80≤CheY-P≤0. 97], within the error of estimation. The optimal adaptation rates which give maximal drift velocities correspond to an average [CheY-P]∼0. 9. In steep gradients, the profile of average CheY-P flattens, and the optimal adaptation rate becomes higher (Figure 6B). The reason why the interval [0. 80≤CheY-P≤0. 97] corresponds to optimal chemotaxis is evident from the profile of motor bias as a function of CheY-P (Figure 6C). The interval [0. 80≤CheY-P≤0. 97] corresponds to the operating range of the motor [0. 95≥mb≥0. 72], where the dependence between mb and CheY-P is approximately linear. In this interval, chemotactic behavior is most efficient in liquid media. The optimal adaptation rate therefore sets the CheY-P level to fit the motor operating range. In steep gradients, the adaptation rate must be high enough to balance the strong excitation and set CheY-P within this optimal interval. In shallow gradients, adaptation must be slow enough to allow excitation, otherwise the cells become adapted before they are able to respond. The effect of varying the [CheR] to [CheB] ratio was studied at fixed [CheB] in three constant-activity gradients N1, N2, and N3 in a liquid medium. The chemotactic efficiency dramatically decreases above [CheR] = 1 (Figure 7), because the resulting higher steady-state CheY-P level produces tumbling behavior. Below [CheR] = 1, chemotactic efficiency decreases slowly for N3, or goes up for the N1 and N2 gradients. The latter effect is caused by a shift of average CheY-P level to the optimal interval, where the chemotactic sensitivity is the highest due to a more optimal fit to the motor operating range. We have further studied the effect of CheB phosphorylation feedback on chemotactic efficiency in a liquid medium. Under the assumption that [CheR] and [CheB] perfectly match each other (A* = 1/3), the CheBp-effect is positive when the adaptation rate is lower than the optimum, and negative when the adaptation rate is higher, in the given gradient (Figure 8A). This effect is caused by the reduction of CheB activity relative to CheR, when the kinase activity A is below the steady-state level (A* = 1/3), as described in the section ‘Model of E. coli Signaling Network’. The average CheY-P level is thus shifted up, which results in a positive or negative effect of CheB phosphorylation, depending on the rate of adaptation (Figure 8B). The positive role of phosphorylation can be significantly increased when the ratio of [CheR] to [CheB] is non-perfect (Figure 8C). For example, 25%-active CheB can significantly counteract the strong negative effect of [CheR] = 1. 25 in the N3 gradient—the drift velocity rises from 1. 8 to 2. 8 µm s−1 (55%). At [CheR] = 0. 75 the effect is not so dramatic, but remains significant—the average drift velocities increase by about 10–15% in all three gradients. This suggests that CheB phosphorylation helps to maintain chemotaxis at fluctuating concentrations of CheR and CheB, when their ratio is not perfect due to gene-expression noise. In the swarm assay in soft agar, bacteria consume an attractant, thereby creating a local gradient, and follow it in the form of a growing ring [54], [55]. We assume that the intensity of the moving gradient remains constant, and use the constant-activity gradient as a simple model for the swarm assay simulation. The constant-activity gradient provides a constant cellular-drift velocity at any distance from the center of the plate. This property allows us to use it as a stationary model of the real moving gradient of attractant. In swarm assays, bacteria move in a labyrinth of agar filaments, with obstacles and traps along the cell' s path. The cell can encounter traps during its run, and stays trapped until it makes the next tumble, as observed by Wolfe and Berg [55]. Therefore, non-adapting cells and non-tumbling mutants form the smallest rings. To simulate motility in such a porous medium as agar, we have introduced a new state of the cell, corresponding to a stop in a trap during a run (Figure 9). The positions of traps are not fixed in space. Instead, it is assumed that each cell encounters traps in an exponentially distributed time series, which mimics the random collisions of the cell with agar filaments. The mean free time between traps is set to 2. 0 s to achieve biologically realistic drift velocities (about 1 µm s−1). While it is trapped, the cell remains stationary until it makes a tumble, whereupon normal run and tumble behavior resumes until the next stop occures [55]. In our model, we assumed that the levels of the adaptation enzymes CheR and CheB vary in a coordinated manner, leaving the [CheR]/[CheB] ratio the same as in the wild type. The ratio of CheR to CheB can be assumed to remain largely fixed because their genes are adjacent and transcriptionally coupled in the meche operon. The adaptation rate in our model is thus proportional to the level of co-expression of CheR and CheB, which will be further denoted as [CheR, CheB]. In order to study chemotactic efficiency at different adaptation rates in agar, we have experimentally measured chemotactic efficiency on swarm plates. In these experiments, CheR and CheB-YFP were co-expressed from one operon under control of a pBAD promoter and native ribosome-binding sites. The pBAD promoter gives expression levels lower or higher than the wild-type value, depending on the strength of arabinose induction. Mean protein levels in the population at a given induction were determined as described in Experimental Methods. Experiment and simulations show that cells with [CheR, CheB] above a certain threshold perform chemotaxis equally efficiently (Figure 10A and 10B). However, the cells with [CheR, CheB] below the threshold have severely impaired chemotactic behavior. According to the simulations, cells with low [CheR, CheB] tend to run without tumbling and stay trapped most of the time. On the other hand, cells with extremely high [CheR, CheB] loose their sensitivity to the gradient and also have poor chemotactic efficiency (Figure S5). This suggests a positive selection for cells with optimal [CheR, CheB] in liquid media—such cells can reach the nutrient source faster and have more available substrates for growth. In contrast, swimming in agar poses mainly negative selection—cells with low [CheR, CheB] are filtered out from the chemotactic population. The limits of motor bias for optimal chemotaxis in agar are also different from those in liquid media. As one can see in Figure 10C, the average CCW motor bias of successful cells is just slightly higher than the steady-state mb0. Cells with higher motor bias would drift faster in liquid media, but not in agar, because the period of time they remain trapped also increases with CCW motor bias. To model swarm assays more realistically, we simulated cell populations with a log-normal distribution of [CheR, CheB] values. The mean (1. 6) and standard deviation (0. 48) are fitted to reproduce the variability of adaptation times observed for wild-type cells [33]: Tad = 311±150 s in response to a 0–10−3 M MeAsp step. The scatter plot of distances travelled by cells along the gradient N2 in a liquid medium shows that a subpopulation with optimal [CheR, CheB] levels drifts more rapidly than other cells (Figure 11A). Simulations in the N3 gradient in agar show that cells with low [CheR, CheB] levels are hindered by agar traps, while other cells drift successfully (Figure 11B). In Figure 11C and 11D the same cells are colored from deep blue to red, according to their [CheR, CheB]. The outer edge of the bacterial ring in a liquid medium contains many blue cells with [CheR, CheB] between 0. 5 and 2. In contrast, the outer edge in the agar contains a uniform mixture of cells with different [CheR, CheB] levels, while deep blue cells with low [CheR, CheB] tend to be left behind. To confirm that chemotactic cells are selected for their [CheR, CheB] levels in swarm plates, cells expressing CheR and CheB-YFP from one operon were taken from two positions in the swarm ring—at the center and at the outer edge—and protein levels in individual cells were determined using fluorescence imaging. The cells collected near the center at a standard agar concentration (0. 27%) have on average lower copy numbers of adaptation enzymes than cells at the outer edge, confirming the predicted selection against low copy numbers (Figure 12A). As expected, in the swarm plates with a reduced agar concentration (0. 20%), the difference between center and outer edge is much smaller (Figure 12B), suggesting that there is no strong selection against low copy numbers in liquid media. It should be noted that agar concentrations below 0. 20% do not produce a stable gel structure, and therefore that is probably the most liquid agar that can be used for swarm plate experiments. Our simulations and additional experiments with a pTrc promoter, which gives much higher basal expression level of [CheR, CheB], show that very high levels of the adaptation enzymes, over 20-fold, can again decrease chemotactic efficiency in agar (Figures S5 and S6).
  In this paper, we present RapidCell—a model of chemotactic E. coli, which allows us to study the effect of chemotaxis network properties on the behavior of large bacterial populations. RapidCell uses a hybrid model for pathway simulation, with mixed algebraic and ODE description instead of a fully stochastic model, AgentCell [19], or a complete system of ordinary differential equations, E. solo [23]. Our model allowed us to dramatically decrease in computational costs. Though many molecular details are skipped or modeled in a rapid-equilibrium (algebraic) approximation, the key steps of the network are reproduced in agreement with up-to-date experimental data. In contrast to detailed single-cell simulation programs which reproduce the noisy behavior of individual cells [19], [56], RapidCell is aimed at predicting the averaged behavior of bacterial populations, and to investigate how it is affected by the signaling network parameters, neglecting the intrinsic noise coming from molecular reactions. However, artificial sources of noise can be further added in the deterministic model of the signaling pathway. In the present version of RapidCell, the noise arises only from rotational diffusion and stochastic switching of the motors. For the receptor cluster simulation, we used the mixed-receptor cluster MWC model [12], [28], [30], which accounts for the observed broad range of sensitivity and reproduces the recent in vivo FRET data [27]. Adaptation is modeled according to the mean-field approximation of the assistance-neighborhood model, with the assumption that the average methylation level of multiple receptors can be represented as a continuous rather than a discrete variable [30]. In contrast to the other reactions, methylation and demethylation are relatively slow and therefore described by an ODE. The ODE is integrated by the first-order Euler scheme to ensure high computational speed of the program, while the time step is chosen as 0. 01 s to keep the simulation error low. Taking into account the available experimental studies on tumble mechanics [31], [57], we use a voting model of run-tumble switching [13], [31], [44]. The model is in a good agreement with experimentally measured run and tumble times. However, more high-resolution experimental data on the interplay among multiple flagella during the run and the tumble would be necessary for a detailed model of run-tumble cellular behavior. There are several types of gradients usually applied in computer models of chemotaxis. The linear gradient arises between stationary source and adsorber, and can often be observed under natural conditions. The Gaussian, another commonly used gradient, appears when a limited amount of molecules is injected into the medium from a micropipette or a similar source [42]. Other gradients that arise from general models of diffusion have hyperbolic or exponential shapes. However, all commonly used gradients have a ‘blind’ zone where receptors are saturated and cells do not respond. When cells drift along these gradients, the average profile of CheY-P changes dramatically, from a steep fall at low concentrations to a weakly stimulated state at high concentrations (Figure 3C). This makes it difficult to compare long-term chemotactic efficiency, because the average CheY-P and drift velocity are non-stable along the gradient. To study chemotaxis systematically, we propose a new—constant-activity—type of gradient. This gradient has the unique property of providing the same CheY-P level and cellular-drift velocity over a wide range of ligand concentrations. The stability of the CheY-P level allows us to study properties of virtual chemotactic cells systematically, and to compare chemotactic behavior over long time periods and concentration ranges. The form of the constant-activity gradient is derived from the MWC model, by formulating the differential equation for the gradient shape which will give a constant rate of receptor free energy change due to ligand binding. In earlier work, the condition of constant chemotactic response was studied using a phenomenological model of ligand binding, with a single dissociation constant KD [48]. The study of Block and co-authors showed that such a model can be simplified, and as a result an exponential ramp of ligand should give a constant response in the range between Cmin = 0. 31KD and Cmax = 3. 2KD, a prediction that was supported by their experiments [48]. In our study, we show that the differential equation for the constant-response gradient proposed in [48] is the result of the MWC model. We further solve this differential equation analytically, and find the exact form of the constant-activity gradient. This gradient grows similarly to the exponential function at moderate ligand concentrations, and increases faster than exponential at low and high concentrations (Figure 2A). Our simulations show that the chemotactic response of the MWC model in the constant-activity gradient remains stable over four orders of ligand concentration—between 0. 1 and 1000KD, in the case when Tsr receptors are fully insensitive to the ligand. However, in the case of (Me) -Asp, the Tsr receptors are able to respond non-specifically to high ligand concentrations, therefore above 100KD the cluster activity drops to zero in a mixed-receptor cluster [12], [27]. However, our simulations of population behavior consider only moderate Asp concentrations, so the cluster activity remains nearly constant in all observed cases. The exponential ramp also gives nearly constant response in the MWC model, but over a much smaller range—between 0. 5 and 3. 0KD, in agreement with [48] and the recent study of Tu et al. [58]. We also show that the apparent dissociation constant KD can be estimated by either the arithmetic or geometric mean of Koff and Kon, but the geometric mean gives a better approximation over a wide range of ligand concentrations. The shape of the constant-activity gradient is also close to a hyperbolic gradient, with the change of variables, KDCx/ (1−Cx) = KD (1/y−1) ∼KD/y, (y = 1−Cx, KD≪1). The hyperbolic gradient arises from simple models of diffusion, when ligand molecules are emitted from a spherical source into the surrounding medium. In nature, such conditions can be observed, for example, in aquatic ecosystems where microalgae leak organic matter attractive for bacteria [59]. This suggests that hyperbolic and exponential gradients with appropriate parameters can be good approximations for the constant-activity gradient. In our model, the adaptation rate is assumed to be proportional to the co-varied concentration of the adaptation enzymes [CheR, CheB], and we use both terms to denote the rate of adaptation. However, increasing expression of the adaptation enzymes may lead to saturation of the adaptation rate at some point, because the enzymes will start working out of saturation kinetics. For these reasons, it is more correct to consider our results in terms of adaptation-rate effects on chemotaxis, whatever the origins of adaptation-rate variability may be. The effect of adaptation rate on chemotaxis agrees in many respects with the results reported in [13] for optimal noise filtering of the chemotaxis signaling system. In their work, the authors demonstrated the existence of an optimal cutoff frequency, an analog of the adaptation rate in our study, for efficient chemotaxis. For a fixed linear gradient, they show the same shape of chemotactic efficiency as a function of cutoff frequency (Figure 3B in [13]) as we found in our simulations (Figure 5A). The authors also show that the optimal cutoff frequency depends on gradient steepness in a linear manner (Figure 5A in [13]), consistent with our results (Figure 5B) for steep gradients. Our simulations in the constant-activity gradient suggest a simple biological mechanism that determines the optimal adaptation rate for a given gradient steepness. Different optimal adaptation rates correspond to a single CheY-P interval, which fits the linear range of the motor-response function. This means that the highest drift velocity in liquid media is observed when the CheY-P level is in the narrow interval fitting the operating range of the motor. In this range, the dependence between CheY-P and mb is approximately linear (Figure 6C). We found that the CheB phosphorylation feedback can have either a positive or negative effect on chemotactic efficiency, depending on how it shifts the average CheY-P level relative to the region of linear motor response. In the case of non-perfect ratio of CheR to CheB, the CheB phosphorylation mechanism can partially counteract the negative effect of unbalanced [CheR]/[CheB], by shifting the average CheY-P towards the optimal region. This confirms that CheB phosphorylation can improve the chemotactic properties of cells with deviations in the ratio of [CheR]/[CheB], as well as in the ratios of other proteins, from the optimum [32]. Chemotactic behavior in liquid media differs from that in agar. We simulated agar effects using traps randomly distributed over time - a cell can encounter traps during its run, and stays trapped until it makes the next tumble, as observed by Wolfe and Berg [55]. This restricts cellular motility—cells that are highly biased towards running remain in traps longer. In agar, the region of optimal motor bias is very narrow and is just above the unstimulated state mb0, because higher bias increases the period of time cells remain in traps. In our model, we did not take into account the growth of a bacterial populations. The typical swarm plate experiments last several hours, and cells grow and divide during the experiment, leading to variations in protein levels and to redistribution of proteins from generation to generation. However, the effect of different adaptation rates in our simulations is clearly visible already within one cell generation over 1000 s of model time (Figure 11B). The selection thus works on a time scale that is shorter than the generation time, which, in our opinion, justifies using a fixed protein distribution. Therefore, the addition of cell growth should not change our results qualitatively. In experiments, daughter cells with sub-optimal levels of CheR and CheB will rapidly fall behind the spreading swarm ring in the vicinity of the division site, while the subpopulation with optimal adaptation rates will be always at the front edge of the ring. In most of our simulations, we assume that the CheR and CheB ratio is constant due to the genetic coupling between the two respective genes, and that cell-to-cell variation in adaptation rates arises from concerted variation in the levels of both enzymes [32]. We also investigated the effects of variation in the [CheR]/[CheB] ratio, which results from translational noise, and affect both the adaptation rate and the steady-state motor bias. In addition to these investigated sources of noise, there is intrinsic noise in the pathway activity which arises from the stochastic nature of (de-) methylation events. The latter sort of noise can also have positive effects on the spreading of cells in a ligand-free medium [56], and even on chemotactic drift in weak gradients [60]. Superposition of variable noise effects on chemotactic efficiency in variable gradients would be an interesting issue for further study. In this work, we have estimated the variability in concerted CheR and CheB concentrations using available experimental data on cell-to-cell variability in adaptation times [33]. We assumed a log-normal distribution for protein concentrations, which also gives a log-normal distribution of adaptation times to a step-wise stimulus from 0 to 10−3 M MeAsp [33]. There are also other experimental estimates of cell-to-cell variation in adaptation times [34] and related simulations [61], but the adaptation rates observed in those experiments were several times higher, presumably due to different culture growth conditions. Our simulations suggest some evolutionary implications. In liquid media with variable food sources and gradient intensities, variability in adaptation times (protein levels) among cells can help the whole population to respond to different gradients more readily, due to positive selection of cells with optimal [CheR, CheB]. In other words, for any given gradient steepness, there will be a subpopulation which has the best [CheR, CheB] to follow this gradient. In contrast, agar poses mainly negative selection on cell populations - cells with low [CheR, CheB] are filtered out from competition, while all other cells travel with approximately equal efficiency. Inspired by the implementation of AgentCell, RapidCell focuses on highly efficient computation of large populations over long periods, keeping cell-response properties consistent with experimental data. The first version of RapidCell allows us to simulate E. coli populations of size 104–105 cells over a time scale of hours, while tracking the signal network dynamics of individual cells. These capabilities permit the modeling of cellular behavior on a macroscopic scale, as in swarm-plate experiments, and the prediction of properties of heterogeneous populations with noisy components of the signaling network.
summary: Chemotaxis plays an important role in bacterial lifestyle, providing bacteria with the ability to actively search for an optimal growth environment. The chemotaxis system is likely to be highly optimized, because on the evolutionary time scale even a modest enhancement of its efficiency can give cells a large competitive advantage. In this study, we use up-to-date experimental and modeling information to construct a new computational model of chemotactic E. coli and implement it in a computationally efficient way to simulate large bacterial populations. Our simulations are performed in a new type of attractant gradient that ensures a constant level of chemotactic excitation at any position. We show that optimal chemotactic movement in a gradient results from a fine balance between excitation and adaptation. As a consequence, steeper gradients require higher optimal rates of adaptation. Simulations demonstrate that the observed intercellular variability of adaptation times, which is caused by gene expression noise, may play a positive role for the bacterial population as a whole, by allowing its subpopulations to be optimally tactic in gradients of different strengths. We further show that optimal chemotactic properties in a porous medium (agar) are different from those in a liquid.
keywords: mathematics
  computational biology/systems biology
  computational biology/signaling networks
title: Dependence of Bacterial Chemotaxis on Gradient Shape and Adaptation Rate