Dependence of Bacterial Chemotaxis on Gradient Shape and Adaptation Rate

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.


Introduction
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 environmentthat is in presence of nutrient gradients-have significant growth advantages, as shown experimentally for different bacterial species [1][2][3][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.
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 10 4 -10 5 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.

Model of E. coli Signaling Network
MWC model. We applied the recently proposed MWC model for a mixed receptor cluster [12,28,29], which accounts for the observed experimental dose-response curves of adapted cells measured by in vivo FRET [27]. An individual receptor homodimer of type r (r = a and s for Tar and Tsr, respectively) is described as a two-state receptor, being either 'on' or 'off'. Receptors form clusters with all receptors in a cluster either 'on' or 'off' together. The clusters are composed of mixtures of Tar and Tsr receptors. At equilibrium, the probability that a cluster will be active is [12]: where F = F on -F off , and F on/off is the free energy of the cluster to be on/off as a whole. For a cluster composed of n a Tar and n s Tsr

Author 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.
receptors, the total free-energy difference is, in the mean-field approximation, F = n a f a (m)+n s f s (m), which is the sum of the individual free-energy differences between the two receptor states where [S] is the ligand concentration, K on=off r is the dissociation constant for the ligand in the on and off state, respectively. The methylation state of the receptor enters via the 'offset energy' e r (m).
The model can also be generalized for binding multiple types of ligand [12,28].
Adaptation model. Adaptation is modeled according to the mean-field approximation of the assistance-neighborhood (AN) model [12,30]. Both CheR and CheB are assumed to bind receptors independent of their activity. A bound CheR (CheB) can (de-)methylate any inactive (active) receptor within the AN. Each bound CheR adds methyl groups at a rate a(12A), and each bound CheB removes methyl groups at a rate bA. Under these assumptions, the kinetics in the AN model are given by We further assume that both enzymes work at saturation: Note that this equation does not imply a first-order reaction mechanism between the adaptation enzymes and receptors-the enzymes work in the zero-order regime. The linear products a(12A)[CheR] (bA[CheB]) mean that a bound CheR (CheB) can only act if the receptor cluster is inactive (active), with probability (12A) and A, respectively.
We further define the relative adaptation rate k: The parameter k indicates the adaptation rate relative to the wildtype adaptation rate V. In the cells with normal steady-state activity (A * = 1/3), the rates and concentrations of the adaptation enzymes are related through b[CheB] = 2a [CheR]. In this work we assume that reaction rates a and b remain unchanged, and the variability in adaptation rate k is caused by variability in [CheR,CheB], provided that they change in a coordinated manner with the fixed ratio: [CheR]/[CheB] = 0.16/0.28 [35]. The latter ODE for methylation is integrated using the Euler method, so that the average methylation level evolves in time as To achieve high computational efficiency in the model, we assumed that the average methylation level m is a continuously changing variable within the interval [0,8], with linear interpolation between the key offset energies: e r (0), 1.0; e r (1), 0.5; e r (2), 0.0; e r (3), 20.3; e r (4), 20.6; e r (5), 20.85; e r (6), 21.1; e r (7), 22.0; e r (8), 23.0, according to [12,30]. Kinase activity. CheA kinase activity is assumed to be equal to the activity of the receptor complex (A). The differential equation for CheY-P is [32] dYp Here  [32,36,37]. The rate of phosphotransfer from active CheA to CheY is much faster than the rate of CheA autophosphorylation (Table S1). Therefore, the concentration of CheY-P is obtained as a function of active CheA from the steady-state equation: In relative units, CheY-P ½ 3 kY ksA kY ksAzkZ Zzc Y , where k s = 0.45 is a scaling coefficient. The relative [CheY-P] = 0,1,3 correspond to fully inactive, adapted and fully active cluster, respectively. The absolute concentration relates to the relative as [CheY-P] abs = 3.1[CheY-P] (mM), following [38].
Time-scale separation. We assume that the rates of ligand binding t l , rates of receptor-cluster conformational changes t k and receptor covalent modification t m are well separated in scales: t l %t k %t m . On our scale (,1 s) the reactions of CheA autophosphorylation, phosphotransfer from CheA to CheY and CheB can be described as a rapid equilibrium state through algebraic equations. The slowest reactions-methylation by CheR and demethylation by CheB-are described through an ODE to reproduce the time scales of seconds and minutes required for adaptation. Table S1 shows the comparative rates of the main reactions.
Model verification. A summary of the parameters used in the model is given in Table 1, and a summary of models and assumptions is shown in Table 2. Along the lines of the MWC model for a mixed receptor cluster [12], we model a cluster of 18 receptors, composed of 6 Tar and 12 Tsr receptors. The catalytic rates a and b were chosen to achieve the proper time scale of adaptation according to in vivo FRET dose-response curves.
As shown previously in [12,29,39], the MWC model for a mixed receptor cluster correctly reproduces the in vivo FRET response amplitudes to step-wise addition and removal of MeAsp [27,40]. We also compare our model output with the published FRET response ( Figure S1A), and show that the simulation is in good agreement with experiment, both for the amplitude and the duration of the chemotactic response. However, the steepness of the adaptation curve after attractant removal can only be roughly described by the existing model of CheB activity, a deficiency which needs to be addressed for more precise modeling in future.
The spatially extended StochSim model gives lower response amplitudes compared to FRET experiments [14]. Comparison of RapidCell and StochSim responses to addition and removal of Asp is shown in Figure S1B. The adaptation rate of StochSim seems very high compared to FRET experiments and RapidCell simulations (k = 8 times higher than the wild-type rate), which suggests that RapidCell will be much more sensitive to gradients than AgentCell [19].
RapidCell also reproduces experimental data on tethered cell stimulation with pulse and step changes of Asp concentration [41] ( Figure S2A and S2B). The adaptation times after a step increase of a-methylaspartate (MeAsp) concentration over three orders of magnitude agree with experimental data reported in [33] ( Figure  S2C).

Model of E. coli Motion
During a run, the cell is assumed to move with a constant speed v = 20 mm/s, while the direction of motion is affected by rotational diffusion [7,42]. After each time step, the running direction is  Table 2. Models used in RapidCell.
Motor switching. The relative concentration of the response regulator [CheY-P] is converted into motor bias using a Hill function [38] (Table 2). Motor bias is the mean fraction of CCW rotation time for a motor: mb = T ccw /(T ccw +T cw ), where T ccw and T cw are the means of exponentially distributed CCW and CW intervals, respectively. The equation gives the frequency of the Poisson process of CCWRCW motor switching. The frequency of reverse switching CWRCCW is l rev = 1/T cw . After each time step Dt, the motor can switch its direction from the present state, according to the current switching frequency l forw(rev) , with probability P forw(rev) = l forw(rev) Dt. Runs and tumbles. Run and tumble events include the complex interplay of filaments in a bundle, the details of which have been investigated experimentally [31,43]. To simulate the run and tumble behavior of a cell with several motors (N = 3-7) we consider the voting model, where the majority of the motors determines the cellular behavior.
Model of voting motors. The cell has N = 5 motors switching independently, and the state of the cell is determined according to a voting model [13,31,44]. The cell switches from 'Run' to 'Tumble', if at least 3 of its 5 motors rotate CW, and from 'Tumble' to 'Run', if at least 3 of the 5 rotate CCW. The choice of N = 5 is arbitrary, and similar results are obtained for N = 3,7 under the condition of majority voting.
For model validation, simulations of cells with N = 3,5,7 motors were carried out ( Table 3). The simulated run times (1.04-1.11 s) agree with the experimental value of 1.2461.16 s [45]. The simulated tumble times (0.26-0.44 s) appear higher than the measured 0.1460.08 s [7,31]. However, the latter study [31] shows that the full tumble time, from bundle breaking in the old run to bundle consolidation in the new is 0.4360.27 s. This estimate of tumble time reflects not only cell reorientation, but also the interplay of flagella and the resulting drop in cell speed, and the voting model reflects specifically this kind of tumble time estimate. The model with 5 motors is used in the following as default.
Tumbling angle. The tumbling angle is distributed according to the probability density function f(H) = 0.5(1+cosH)sinH, 0,H,p [46,47], with M(H) = 67.5u which is close to the experimental measurement of 68u [7], and the corresponding shape of the function ( Figure S3).

Model of the Environment
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, x max )6 (0, y max ), and the circular domain within (0, r max ), with x max = y max = 2r max = 20 mm.
The constant-activity gradient. The gradients used in chemotaxis modeling are usually linear, Gaussian or exponential [19,23]. However, in these gradients the signal is non-constant, which means it is strong at low attractant concentrations, and weak at high concentrations due to receptors saturation. Such a non-uniform distribution of the signal makes it difficult to estimate chemotactic efficiency over long time intervals-cells soon become 'blind' because receptors are saturated, and chemotactic drift decreases.
According to the MWC model, an increase in ligand concentration DS causes an initial rise in receptor free-energy difference Using the Taylor-series approximation, leads us to the following approximation for free energy per concentration change: Simplified solution. The denominator in Eqn. 13 can be simplified by assuming and the unknown K * can be found from equation which gives two alternative estimates for K * :K Ã~K on zK off 2 and K Ã~ffi ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi K on K off p , i.e. the arithmetic and geometric means of K on and K off .
At zero or relatively low ligand concentrations, the geometric mean has a high impact in Eqn. 17, and is preferable as an estimate. Indeed, in earlier work it was earlier referred to as the apparent dissociation constant K D of ligand binding [14]. However, at high concentrations, the arithmetic mean will have a higher impact in Eqn. 17, so it can be used as an alternative estimate. Our simulations indicate that within four orders of aspartate concentration the geometric mean serves as the best estimate of K * ( Figure S4). Taken together, the energy difference is approximated by describes the unknown function S(x), which will give the 'constantactivity' gradient shape. The function S(x) will give a constant change of energy difference C per length unit dx of cellular path along the gradient. In other words, such a shape of gradient will give a constant cluster activity at any ligand concentration.
Within the accuracy of a constant term, the latter differential equation was previously used by Block and Berg in [48], who derived it assuming that receptor occupancy is proportional to S/ (S+K D ), with a single K D for active and inactive receptors. The authors assumed that the chemotactic response is proportional to the change in receptor occupancy [27,48]. They simplified this equation to reduce the variability of the 1 KDzS ð Þ 2 term, leading to the exponential form of the solution.
However, we can solve Eqn. 18 analytically: Þ is the constant of integration, determined by the initial condition S(0). The condition S(0) = 0 gives the following chemoattractant function: Constant-activity gradient of Asp. In the case of aspartate (K on = 12, K off = 1.7, K * = 4.52 mM), the S(x) function reads: Our simulations demonstrate that this form of constant-activity Asp gradient gives a constant cluster-activity response with reasonably good precision (see Results). Gradient steepness. A cell swimming with speed v along the axis X from the point (x = 0) senses the monotonically increasing function S(x) and a constant change in receptor free energy dE=dt~Cdx=dt~Cv ð22Þ per second, which is defined as the steepness of the constant-activity gradient.
Limiting condition. Note the necessary condition for Eqn. 20 to avoid singularity and negative concentrations. It sets the upper limit CvC max~K xmax for the gradient steepness C within the domain (0, x max ). For example, within a domain of size x max = 10 mm, the maximum steepness of a gradient of aspartate is Cv = 2.28/x max v = 4.56610 23 .
Constant-activity and exponential ramps. In contrast to spatial gradients, which direct the cellular motility in a certain direction, time ramps are used to study the chemotactic response of tethered cells [41,48].
The constant-activity ramp of Asp was simulated according to Eqn. 20: with simulation time T max = 1000 seconds. The resulting value of C gives the maximum ligand concentration S(T max ) = 9999K * .
The exponential ramp was simulated as: after 200 s of adaptation to the initial stimulus 0.31K D , following the model and experiments of [48]. The concentration profiles are shown in Figure 2A.
Constant-activity gradient simulations. The constantactivity gradient (Eqn. 20) has an intensity C~K xmax , and the domain has a rectangular (0, x max )6(0, y max ) or circular (0, r max ) shape. The gradient has its minimum S = 0 at x = 0 (or r = 0) and reaches its maximum S = 999K * at x = x max (or r = r max ) ( Figure 3A). In most simulations we used the circular gradient S(r), and the cells start swimming in random directions from the center r = 0.
Chemotactic efficiency. Chemotactic efficiency was estimated as the average drift velocity of a cell population, measured between 200 and 500 s of simulation time, in the three basic constant-activity gradients N1, N2, N3. As shown in Figure 4, within this interval the average CheY-P level of cells is constant, and the drift velocity can be accurately measured by a linear fit.
Population behavior. The population behavior in the absence of attractant fits the diffusion equation AEr 2 ae = 4Dt. Simulations give a diffusion coefficient D = 2.56610 26 cm 2 s 21 , in agreement with the experimental D = 2.5-3.8610 26 cm 2 s 21 (see [45] and the review of other published values therein).

Program RapidCell
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 Dt = 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). RapidCell is platform-independent and runs as a console application. Its implementation provides a computational speedup of 8000 times compared to AgentCell (based on StochSim without receptor coupling), and approximately 100 times compared to BCT. It enables simulations of up to 100,000 cells to be completed within a time frame of hours using a desktop computer with comparable CPU power and RAM to those mentioned above.

Experimental Methods
Strains and plasmids. E.coli strain RP2867 (tap cheR cheB) is a derivative of RP437 [51]. Plasmid pVS571 encodes cheR and cheB-eyfp as parts of one operon under control of a pBAD promoter and native ribosome binding sites. The insert cheR cheB-eyfp was recloned with SacI and XbaI from the plasmid pVS145 which was constructed by cloning a PCR-amplified fragment containing cheR upstream of cheB-eyfp in the pVS138 plasmid [52] using a SacI site introduced by the upstream PCR primer and a HindIII site in cheB.
Swarm experiments in soft agar plates. Tryptone-broth (TB; 1% tryptone, 0.5% NaCl) soft agar plates were prepared by supplementing TB with 0.27% agar (Applichem), 34 mg ml 21 chloramphenicol, and indicated concentrations of arabinose. Cells  were inoculated from fresh colonies grown on LB agar plates. Swarm assays were performed at 34uC for 10 hours or at 30uC for 17 hours. Following incubation, photographs of plates were taken using a Canon EOS 300 D camera, and subsequently analyzed with ImageJ (Wayne Rasband, NIH) to determine the diameter of the swarm rings.
Quantification of gene expression. For quantification of mean expression levels of the fluorescent reporter protein CheB-YFP, cells were grown in liquid TB medium supplemented with 34 mg ml 21 chloramphenicol, and indicated concentrations of arabinose. Fluorescence was determined in a population of cells using flow cytometry on a FACScan (BD Biosciences) equipped with a 488 nm argon laser [32,52]. The autofluorescence background was measured for control cells and subtracted from all values. Single-cell levels of fluorescent reporter proteins in swarm assays were measured by fluorescence imaging on a Zeiss AxioImager microscope and quantified with an automated custom-written ImageJ plugin [52].
To calibrate the fluorescence intensity in FACS and imaging data, a PerkinElmer LS55 luminescence spectrometer was used to determine the absolute number of reporter proteins in control cells. The cells were sonicated with a Branson Sonifier 450 until complete lysis was achieved and YFP fluorescence was measured at 510 nm excitation and 560 nm emission. Sonicated cells without a fluorescence reporter were used as a negative control, and their autofluorescence was subtracted from all values as background. A solution of purified YFP of known concentration, determined by Bradford assay and absorbance measurement by a Specord205 spectrophotometer (Analytik Jena), was used to produce a calibration curve, relating fluorescence to molecule number. Cell number in 1 ml culture was counted using a Neubauer counting chamber, and cell volume was determined by measuring cell width and length by imaging. These values from one culture were used to provide a conversion factor from FACS or imaging values to single-cell protein levels.

Chemotaxis in Different Gradients
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.
Response of the MWC model to ramps. It was previously shown that tethered cells respond with constant strength to an exponentially rising gradient of MeAsp, in the range between 0.31 and 3.2K D [48]. We simulated the response of the MWC model to increasing ramps of Asp in the exponential and constant-activity form (Figure 2A). Indeed, the exponential ramp gives nearly constant response between 0.5 and 3.0K * , consistent with the model of [48].
However, the constant-activity ramp results in a chemotactic response that remains approximately constant over three orders of ligand concentration-between 0.1 and 100K D ( Figure 2B). If Tsr is non-sensitive to the ligand, constant activity remains up to 1000K D . However, since Tsr receptors are able to respond to aspartate non-specifically, the activity drops to zero, as previously shown for a mixed-receptor cluster [12,27].
Chemotactic efficiency of cell populations in different gradients. To study chemotactic efficiency in common gradients that arise from general diffusion models, we simulated chemotactic motility in linear and Gaussian gradients ( Figure 3A), and compared them with the constant-activity gradient. The chemotactic efficiency was estimated by the average drift velocities of populations consisting of 1000 identical cells. In Figure 3B, one can see that in the linear and Gaussian gradients the drift velocity decays after about 400 and 800 s, respectively, indicating that cells loose sensitivity due to receptor saturation. In contrast, the constant-activity gradient keeps the drift velocity constant at any point ( Figure 3B), as expected.
This population behavior can be explained by the intracellular CheY-P levels of the cells in these gradients. Gaussian and linear gradients result in a strong excitation at low attractant concentrations, and poor excitation at high concentrations ( Figure 3B). In contrast, the constant-activity gradient produces an approximately constant level of CheY phosphorylation across the cell population ( Figure 3C). These two unique properties of the constant-activity gradient-constant drift velocity and constant average CheY-Pfavor this gradient as a reliable in silico assay to study the chemotactic motility of cells.
Average CheY-P level in the constant-activity gradients. Simulation of cell populations in the constantactivity gradients N1, N2 and N3 demonstrate that the average CheY-P level depends on gradient steepness and remains stable Relative adaptation rate k = 1. The cell parameters are as described in Table 1. doi:10.1371/journal.pcbi.1000242.g004 over long time intervals (Figure 4). These three gradients were used further, as default, to measure chemotactic efficiency under different test conditions.

Optimal Adaptation Rates in a Liquid Medium
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 (nonadapting 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.88610 23 , ( 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 mm s 21 ) [53]. Our simulations indicate that experimental celldrift velocities are inlikely to exceed 15 mm s 21 , corresponding to an extremely steep and short-scale gradient. In very weak gradients, the drift velocity can be as low as 2.5 mm s 21 , still distinguishable from the non-chemotactic cellular drift (0.8 mm s 21 ). 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.

Effect of [CheR] to [CheB] Ratio on Chemotactic Efficiency
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.

Effect of CheB Phosphorylation on Chemotactic Efficiency
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 mm s 21 (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 geneexpression noise.

Swarm Plate Simulations
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 constantactivity 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, nonadapting 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 Gray bands indicate the optimal CheY-P interval and the corresponding operating range of the motor. The cell parameters are as described in Table 1 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 mm s 21 ). 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].

Optimal [CheR,CheB] in Agar-Experiments and Simulations
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 mb 0 . 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.

Swimming in a Liquid Medium and Agar with a Log-Normal Distribution of [CheR,CheB]
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]: T ad = 3116150 s in response to a 0-10 23 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.

Measurement of [CheR,CheB] in Individual Cells in Different Parts of Swarm Rings
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 Kinase-dependent CheB activity means that CheB works more weakly at A,1/3, and thus the average [CheY-P] is higher than the level obtained for constantly active CheB. Such a shift improves chemotaxis at low adaptation rates, but reduces it at high rates. The optimal range of CheY-P is shown by the gray band. (C) Drift velocities at variable [CheR] and variable CheB activity and fixed [CheB] (0.28 mM, wild type). Solid, dashed and finely dashed lines indicate 100%, 50% and 25% active CheB, respectively. Adaptation rate k = 1, other cell parameters as described in Table 1. doi:10.1371/journal.pcbi.1000242.g008 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   Figure 4 of [55]. The cross shows the drift velocity of non-adapting cells, from Figure 6 of [55]. The cell parameters are as described in Table 1 [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).

Discussion
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 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 newconstant-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 K D [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 C min = 0.31K D and C max = 3.2K D , 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 1000K D , 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 100K D 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.0K D , in agreement with [48] and the recent study of Tu et al. [58].
We also show that the apparent dissociation constant K D can be estimated by either the arithmetic or geometric mean of K off and K on , 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, K D Cx/ (12Cx) = K D (1/y21),K D /y, (y = 12Cx, K D %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 mb 0 , 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 23 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 10 4 -10 5 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. Figure S1 Comparison of the RapidCell network response with experimental and simulated data. (A) FRET experiment and RapidCell simulation of cell response to a step-wise stimulus of MeAsp. The initial ambient concentration is zero; at t = 80 s 30 mM MeAsp is added and removed at 480 s. The best fit by RapidCell is obtained with an adaptation rate of k = 0.5, corresponding to the temperature T = 20uC at which the FRET experiments were carried out. At T = 30uC, the fitted adaptation rate will be k = 1.0 (V.Sourjik, unpublished data). (B) StochSim and RapidCell simulations of cell response to a step-wise stimulus of Asp. The initial ambient concentration is zero; at t = 20 s 3.5 mM Asp is added and removed at 70 s. The best fit by RapidCell is obtained with an adaptation rate of k = 8 -a very rapid rate of adaptation. The StochSim simulations were carried out with a coupled model (Shimizu et. al, 2003), consisting of 65665 square receptor lattice with coupling energy E J = 23.1 kT. Found at: doi:10.1371/journal.pcbi.1000242.s001 (0.30 MB TIF) Figure S2 Comparison of the RapidCell network response with experimental data on tethered cells. (A) Simulation of CCW motor bias response to a short pulse of attractant. The initial ambient concentration is zero; at t = 5 s 1.0 mM Asp is added for a 0.35 s interval; solid line -simulations (the best fit is obtained with an adaptation rate of 2.0), circles -experimental data (Segall et. al., 1986). (B) Simulation of CCW motor bias response to a step-wise stimulus. The initial ambient concentration is zero; at t = 1 s 0.075 mM Asp is added; solid line -simulations, circlesexperimental data (Segall et. al., 1986). The best fit is obtained with an adaptation rate of 5.0. (C) Adaptation times to a step increase of MeAsp from zero ambient level, obtained in simulations (solid line) and in experiments (Berg and Tedesco, 1975) (circles). In the simulations, the dissociation constants used were K a off = 0.02 mM and K a on = 0.5 mM (Keymer et. al., 2006). The best fit is obtained with an adaptation rate of 1.3.