Sequential modification of bacterial chemoreceptors is key for achieving both accurate adaptation and high gain

Many regulatory and signaling proteins have multiple modification sites. In bacterial chemotaxis, each chemoreceptor has multiple methylation sites that are responsible for adaptation. However, whether the ordering of the multisite methylation process affects adaptation remains unclear. Furthermore, the benefit of having multiple modification sites is also unclear. Here, we show that sequentially ordered methylation/demethylation is critical for perfect adaptation; adaptation accuracy decreases as randomness in the multisite methylation process increases. A tradeoff between adaptation accuracy and response gain is discovered. We find that this accuracy-gain tradeoff is lifted significantly by having more methylation sites, but only when the multisite modification process is sequential. Our study suggests that having multiple modification sites and a sequential modification process constitute a general strategy to achieve both accurate adaptation and high response gain simultaneously. Our theory agrees with existing data and predictions are made to help identify the molecular mechanism underlying ordered covalent modifications.

M ost post-translational regulatory processes involve reversible covalent modifications (phosphorylation/ dephosphorylation, methylation/demethylation, etc.) of key proteins catalyzed by enzymes (kinase/phosphatase, methyltranferase/methylesterase, etc.). Instead of having only a single site of modification, many regulatory proteins such as histones, p53, RNA polymerase II, tubulin, etc. have multiple modification sites 1 . The multiple modification sites allow a single regulatory protein to have complex functions depending on combinations of different modification processes 2 . For example, the histone proteins have multiple covalent modification sites of different types (methylation, acetylation, phsophorylation, etc.) and the different combinations of the multiple modification sites are thought to code for different gene expression patterns in different cells 3 . However, how this combinatorial molecular code works, i.e., how it is encoded and decoded, remains poorly understood 4 .
One well-studied multisite regulatory protein is the cyclindependent kinase inhibitor Sic1 in Saccharomyces cerevisiae (yeast). Sic1 has more than six phosphorylation sites whose main function is regulating the timing of the G1/S transition in yeast cell cycle 5,6 . Huang and Ferrell 7 first suggested that the response sensitivity can be enhanced by having multiple modification sites. However, Gunawardena 8 pointed out that other effects such as substantial disparities in enzyme efficiency among different sites are also needed in making a sharp switch. Later work by Salazar and Hofer 9 showed that a random phosphorylation process among the different sites gives rise to a shallow but rapid response while sequential processing gives rise to a steeper but slower response. Though much progress have been made, dynamics and functions of multisite modification in Sic1 remain not fully understood.
In this paper, we focus on a relatively simpler signaling system, bacterial chemotaxis 10 , where multisite modification has an important role in adaptation 11 . Adaptation is an important general biological behavior that allows a living system to adjust its internal state in response to changes in its environment so that it can return to a set activity level after a fast response to a persistent change in the external stimulus 12 . In bacterial chemotaxis, a chemoreceptor has multiple methylation sites. The kinase activity of a chemoreceptor is determined by chemoeffector ligand concentration (external stimulus) as well as the receptor methylation level (internal state)-a higher attractant concentration leads to a lower kinase activity and a higher methylation level leads to a higher kinase activity. Adaptation in bacterial chemotaxis is achieved by a feedback mechanism in which the receptor methylation level (internal state of the receptor) is controlled by a methyltransferase CheR and a methylesterase CheB that act at a time scale much longer than the response time to a change in external stimulus (ligand concentration). The catalytic efficiencies of CheR and CheB depend on the receptor activity, which form the feedback mechanism for adaptation [13][14][15][16] . However, despite the general consensus on the importance of a negative feedback mechanism for accurate adaptation in bacterial chemotaxis, the detailed receptor methylation/demethylation kinetics among the multiple methylation sites remain unclear.
How do different methylation kinetics, random or sequential, affect adaptation accuracy and response gain? What are the benefits to have multiple modification sites? In this paper, we address these questions by systematically investigating effects of different multisite modification processes, from purely random to strictly sequential, on system-level functions such as adaptation accuracy and signal amplification. Our theoretical findings allow us to infer the multisite modification dynamics from existing experimental data. More importantly, our study leads to specific suggestions of future experiments to determine the molecular mechanism controlling the multisite modification dynamics.

Results
Modeling multisite modification dynamics. In previous modeling studies of the receptor methylation (demethylation) reactions, the microscopic methylation state of the receptor (μ) has been ignored. Here, we consider the transitions between the 2 M = 16 (M = 4 is the total number of modification sites) individual microscopic methylation states of a receptor explicitly. As shown in Fig. 1, all the states μ are grouped (column-wise) by their total methylation level m ¼ P 4 j¼1 μ j , so a methylation (demethylation) reaction moves the current state to another state in the column to the right (left). However, among the multiple states in the next column which one does it transition to? And at what rate? Here, we consider two cases, one special and one general, as shown in Fig. 1a, b, respectively.
For the special case of strictly sequential modification, which is implicitly assumed in previous models 17,18 , the methylation and demethylation processes follow the same sequence (in opposite directions) among the five states shown in Fig. 1a. Following previous work [13][14][15] , the negative feedback control is implemented by only allowing methylation (demethylation) of the inactive (active) receptors respectively. If we define k þ m and k À m as the average methylation and demethylation rates for all receptors with the same total methylation level m, this negative feedback mechanism leads to: where a h i m is the average activity of receptors with methylation level m and the kinetic rates, k + and k − , are proportional to CheR and CheB concentrations, respectively.
In the general case, when site j − 1 is methylated (μ j−1 = 1), the methylation rate for the next site in the sequence j in state-μ is given by the same sequential methylation rate k þ m as in the strictly sequential case described above. However, when site j − 1 is not methylated (μ j−1 = 0), methylation of site j can still occur via the random methylation process albeit with a smaller rate k þR where η is a parameter (0 ≤ η ≤ 1) characterizing the randomness in the methylation process. Combining these two possibilities, the site-specific methylation ratek þ j for site j can be written as: Similarly, demethylation of site j depends on whether site j + 1 is demethylated, and the site-specific demethylation ratek À j for site j can be written as: where k ± m in Eqs. (2) and (3) are the sequential methylation and demethylation rates given by Eq. (1). To describe modifications of the boundary states, i.e., the fully unmethylated state (m = 0) and the fully methylated state (m = M), we introduce a forward initiator with μ 0 = 1 for methylation of the j = 1 site in Eq. (2) and a reverse initiator with μ M+1 = 0 for demethylation of the j = M site in Eq. (3).
Despite the same feedback mechanism given by Eq. (1), dynamics of receptor methylation level depends on the degree of randomness (η) in the multisite modification process. In this paper, we investigate consequences of different multisite modification schemes, from sequential to random, by studying behaviors of the standard model of bacterial chemotaxis for different values of 0 ≤ η ≤ 1. In particular, we study how the adaptation error ξ and the response gain Γ are affected by η. Details of the full standard model framework for studying bacterial chemotaxis and the precise definition of ξ and Γ are given in the Methods section.
Sequential modification is essential for perfect adaptation. In general, the adapted activity level a h i A ð½LÞ is a function of the ligand concentration [L]. Adaptation is deemed perfect if a h i A is a constant independent of [L]. For a general case 0 < η < 1, we can determine the adapted activity by solving the full model numerically using Monte-Carlo method (see Supplementary Methods for details). However, simple analytical equations for 〈m〉 can be found for the extreme cases η = 0 and η = 1, which provide insight on a key condition for perfect adaptation: where 〈a〉 is the average receptor activity and the term ϵ comes from the boundary effects at m = 0 and m = M (see Supplementary Note 1 for details of the derivation). For the case of purely sequential methylation (η = 0), the right hand side of Eq. (4) has the remarkable property of only explicitly depending on 〈a〉 but not on 〈m〉 or [L]. As a result, the adapted activity a h i A % k þ =ðk þ þ k À Þ is independent of [L], i.e., perfect adaptation 17,18 . The independence of the methylation rate dhmi dt on m is that only one modification site is available for methylation or demethylation per receptor at any given time when modification reactions are sequential. Figure 2a illustrates the adaptation process in response to a series of step increase in ligand concentration The solid line represents the adapted activity obtained by setting the right hand side of Eq. (4) to zero. The dashed curves represent the activity as a function of 〈m〉 for different values of [L]. Upon a sudden increase of [L], say from [L] 1 to [L] 2 , the system first responds by decreasing its activity as represented by the downward arrow (blue) as illustrated in Fig. 2a. This altered activity triggers the adaptation mechanism that slowly increases m, causing the system to follow the upward arrow (green) along the dashed line for [L] 2 until it reaches the adapted activity level that is roughly independent of [L]. The fundamental reason for perfect adaptation is that a h i A is independent of 〈m〉, i.e., the solid line in Fig. 2a is flat for a large range of 〈m〉.
For the case of random methylation (η = 1), all the available modification sites are equally accessible. Therefore, the methylation and demethylation rates are proportional to (M − m) and m respectively as given in Eq. (5). As a result, the adapted activity has a simple linear dependence on adapted methylation a h i A ¼ ðM À m h i A Þ=M as shown in Fig. 2b. An increase of ligand concentration from [L] 1 to [L] 2 triggers an immediate response (a drop in activity) followed by a slow adaptation process that leads the system to a different adapted activity level. The inaccurate adaptation for η = 1 is caused by the explicit dependence of hai A on hmi A , i.e., the solid line in Fig. 2b is tilted. It is easy to see that the dependence of hai A on hmi A occurs for all η ≠ 0.
Results from direct Monte-Carlo (MC) simulations of 〈a〉 subject to a series of step increases in concentration [L] (×10 fold increase for each step) as shown in Fig. 2c for sequential and Fig. 2d for random methylation schemes support our analysis shown in Fig. 2a, b, respectively.
For sequential modification (η = 0), the adaptation error ξ (see Methods section for its definition) is proportional to the probability of receptors in the extreme (boundary) methylation states m = M & 0 and we have: where c 1 and b are constants, ξ 0 is the error from the m = 0 state, and α (<0) is the free energy change for adding a methyl group to the receptor (see Eq. (10) in Methods section for the definition of α). Equation (6) shows that ξ decreases exponentially with |α| before saturating to ξ 0 . For random modification (η = 1), the adaptation error has contributions from the whole range of methylation levels 0 ≤ m ≤ M and we have: which only decreases with |α| algebraically (see Supplementary Note 1 for details of the derivation for Eqs. (6) and (7)). The different dependence on |α| given in Eqs. (6) and (7) are verified by direct simulations (Supplementary Fig. 4a), which clearly show that sequential modification reduces adaptation error much more efficiently than random modification.
The tradeoff between response gain and adaptation accuracy. Besides the adaptation error ξ or equivalently the adaptation accuracy ξ −1 , another important property of the system is its response gain Γ, which measures the sensitivity of the system in response to a change in external signal (see Methods section for the definition of Γ).
As shown in Eqs. (6) and (7), adaptation accuracy can be increased by increasing |α|, but what happens to the gain Γ? Interestingly, increasing |α| leads to a reduced gain independent of whether the modification dynamics is sequential or random ( Supplementary Fig. 4b). The reason is that for a larger value of |α|, individual receptors in the receptor cluster in the adapted state will have activities further away from the adapted mean value 〈a〉~1/2-either closer to 0 or closer to 1-where the sensitivity (gain) is lower (see Supplementary Methods for details). It is worth noting that this dependence of Γ on |α| is due to the discrete methylation level of individual receptor, which is only captured by the Ising model [19][20][21] but not in the simplified Monod-Wyman-Changeux (MWC) model [22][23][24] .
The tradeoff or anti-correlation between response gain Γ and adaptation accuracy ξ −1 is a general property of the signaling pathway. Besides the extreme cases (η = 0 and η = 1) considered so far, this tradeoff between Γ and ξ −1 exists for all intermediate cases of methylation dynamics with 0 < η < 1. As shown in Fig. 3a, the gain is almost unaffected when we change the value of η while keeping the other parameters constant, but the corresponding adaptation accuracy ξ −1 decreases with η. On the other hand, when we tune other parameters to maintain a high accuracy (e.g., by increasing |α|), the corresponding gain goes down with η as shown in Fig. 3b. Therefore, for a more random methylation scheme (a larger value of η), the tradeoff between ξ −1 and Γ means that one is enhanced at the expense of the other.
Accurate adaptation and high response gain represent two of the most desirable but opposing properties of biological signaling systems, i.e., to resist changes in the environment by adaptation and to respond to weak signals. This accuracy-gain tradeoff is related to the fluctuation-dissipation relationship established in equilibrium systems 25 . Next, we show how this tradeoff can be , e.g., from [L] 1 to [L] 2 , the system responds quickly by decreasing its activity (blue arrow) from the old adapted state (solid red circle) to the maximum response state (hollow red circle) without changing 〈m〉. This initial response is followed by the slow adaptation dynamics (green arrow) along the dashed line until the new adapted state is reached. Direct Monte-Carlo simulation results of the average activity 〈a〉 in response to a series of step increase in methyl aspartate concentration over 7 orders of magnitudes are shown for c Sequential (η = 0) and d Random (η = 1) modification processes. The step changes in stimulus is shown in Supplementary Fig. 8. The sequential modification process leads to a much higher adaptation accuracy than the random modification process.  Fig. 3 The tradeoff between the adaptation accuracy and response gain. a As η increases, the adaptation accuracy ξ −1 , red squares, decreases while the signaling gain Γ, green circles, remains roughly constant. b When parameters are tuned to keep the accuracy ξ −1 roughly constant for different values of η, the corresponding gain Γ decreases with η. The range of stimulus is set by ½L min ¼ 1 μM and ½L max ¼ 100 mM.
attenuated by having multiple modification sites and sequential modification.
Sequential modification attenuates the accuracy-gain tradeoff.
Why are there multiple modification sites in a regulatory protein or a receptor? How is the performance of the system enhanced by having multiple modification sites? Here, we investigate how having multiple modification sites affects the gain Γ and accuracy ξ −1 for different modification dynamics (sequential versus random). For the sequential modification dynamics, the adaptation error comes from the receptor populations with the extreme (boundary) methylation levels m = 0 or m = M. As the probability P M of reaching the boundary state m = M decreases exponentially with M, the adaptation error in the sequential modification model (η = 0) should decrease strongly (exponentially) with M as given in Eq. (6). For the random modification dynamics, the adaptation error comes from all methylation levels and the reduction of adaptation error with increasing M is much smaller~1/M as given in Eq. (7).
We studied the dependence of the performance of the system on M systematically by computing Γ and ξ for a random set of parameters for M = 1, 2, 3, 4 in our models with η = 0 and η = 1. The results, as shown in Fig. 4, clearly demonstrate the general accuracy-gain tradeoff, i.e., the inverse dependence of Γ and ξ −1 in all cases studied. However, there are significant differences between the sequential and random modification cases. For sequential modification (η = 0), the tradeoff curve is lifted significantly as M is increased as shown in Fig. 4a. In fact, the threshold lines (solid lines in Fig. 4), which are just fits to the highest performing points for each value of M, follow an approximate form: where C 0 (M) measures the overall performance of the system with sequential modification (η = 0). As shown in the inset of Fig. 4a, C 0 (M) increases significantly (linearly) with M. In contrast, as shown in Fig. 4b, the threshold lines in the random modification case follow a much more gradual curve: where the overall performance C 1 (M) for the random modification system (η = 1) has only a weak dependence on M (see inset in Fig. 4b).
The significantly different dependence of the accuracy-gain tradeoff relationship on M for η = 0 and η = 1 clearly shows that having multiple modification sites can ease the accuracy-gain tradeoff in general but the effect is significant only when the modification dynamics are sequential.
Comparisons with existing experiments. In this section, we discuss specific model results that can be directly compared with existing experiments. The E. coli chemoreceptor Tar has four methylation sites at residues 295, 302, 309, and 491, which are labeled by numbers 1-4, respectively. Protein methylation in eukaryotic cells is usually associate with lysine or arginine residues. However, glutamate is the most common residue for methylation in E. coli 26 , and bacterial chemotaxis receptors in general are methylated in glutamate residues; or in glutamine residues that were posttranslationally deamidate to glutamates by CheB. Residues 1-3 are seven residues apart from each other, along the same α helix, whereas residue 4 is located on another helix 27 as illustrated in Fig. 5.
Experimental results 28 indicate that methylation of sites 1, 2, and 3 depends on each other in reverse order, i.e., site 3 is methylated first, followed by site 2 and then site 1, and that residues 316 and 498 affect the methylation of site 3 and 4, respectively. Structural models 29 of the receptor modification are consistent with the methylation rate depending on a residue seven residues away in the C-terminal direction. The initiator residues,  ig. 5 Illustration and notation for the Tar receptor. The methylation state of site i (=1, 2, 3, 4), green circles, is described by a binary variableμ i (1methylated; 0-unmethylated). The sequential methylation and demethylation processes among sites 1-2-3 are shown by the red and orange arrows. The two initiator sites (316 and 288), blue circles, are described by two binary numbersh The receptor methylation state is described by six binary numbers,μ 4 ðh þ 3 Þμ 3μ2μ1 ðh À 1 Þ. In our notation, a methylation site is modifiable when it is labeled by x and a specific value (0 or 1) is assigned when it is fixed by mutation. The two initiator residues are given byh þ 3 andh À 1 -whenh þ 3 ¼ 1, methylation at site 3 (μ 3 ) becomes enhanced; whenh À 1 ¼ 0, demethylation of site 1 (μ 1 ) becomes enhanced. Otherwise whenh þ 3 ¼ 0 orh À 1 ¼ 1, the initial methylation of site 3 or the initial demethylation of site 1 are controlled by the slow random methylation or demethylation processes.
We note that though the existence of theh þ 3 site is supported by experiments [28][29][30] , the initiator siteh À 1 for demethylation is introduced here hypothetically according to the close relationship between the two enzymes CheR and CheB as suggested in a study by Djordjevic et al. 31 , which stated that "structural similarity between the two companion receptor modification enzymes, CheB and CheR, suggests an evolutionary and/or functional relationship" and "the proposed receptor interaction clefts occur on different faces of the β-sheet in CheB and CheR. Topological differences in the structures of CheR and CheB may be reflective of their functionally antagonistic interactions with the receptors".
We first study the adaptation accuracy and response gain from our model and compare them with available experiments. For wild-type (wt) cells with both CheR and CheB, though there is no direct measurements of the methylation dynamics, there have been detailed experimental studies of the in vivo kinase activity dynamics in response to a wide range of stimuli [32][33][34] , which can be compared with our model to determine the response gain Γ and adaptation accuracy ξ −1 .
In ref. 32 the relative sensitivity, S r , is defined as the fractional change in the FRET signal divided by fractional change in stimulus S r ¼ ΔFret=Fret Δ½L=L $ g, where the FRET signal is proportional to the kinase activity and g is the integrand of Eq. (17). From the experimental data on S r (first peak in Fig. 3 of ref. 32 ) and our model, we can estimate the value for the gain Γ for Tar: 4.5 ≲ Γ ≲ 5.
From the measured adapted activity for different background ligand concentrations as plotted in Fig. 1B in the paper by Neumann et al. 33 , we obtained the value of adaptation accuracy for Tar in response to methyl aspartate with a maximum concentration ½L max = 5-10 mM to be roughly in the range: 2.3 ≤ ξ −1 ≤ 3.5.
These estimated values of gain and adaptation accuracy, shown as the black diamond in Fig. 4a, suggest that methylation dynamics should be mostly sequential. Quantitatively, from our model and by using the values of ξ and Γ, we can determine the range of the effective randomness parameter for Tar: 0.05 ≤ η ≤ 0.13 (see Supplementary Methods and Supplementary Fig. 2 for more details on comparison between simulation results and experiments).
Next, we study the methylation profiles in CheB − mutants by using our model and compare them with existing experiments. Different modification dynamics, random or sequential, lead to qualitatively distinct mean methylation profiles at a given time. Denote p j (t) the probability of site j being methylated at time t. For a purely random modification dynamics, all p j (t) should be the same. However, for methylation dynamics that are dominated by sequential modification, the methylation pattern among different sites follows certain distinctive patterns.
The methylation dynamics among different sites in the demethylase CheB − mutants were studied by the Koshland lab more than 20 years ago 28,30 , by mutating the residuesμ 4 ,h þ 3 ,μ 3 , μ 2 in Fig. 5, and the initiatorh þ 4 . After being methylated with tritiated SAM, the receptors were cleaved and the extent of methylation of each site was determined by high-performance liquid chromatography. The methylation rates were calculated in arbitrary units, reproduced here in Table 1. As the absolute values are not available, we can only analyze the relative methylation ratios of the different sites and mutants.
It is useful to compare simulations with k − = 0 with the CheB − mutants to isolate the effects of sequential methylation. Specifically, these mutants are, besides the wild-type receptor (EEQE), the mutant receptors EEDE and EEEE. In these strains, site 3 (E309) can be methylated when occupied by Glutamate (E) residues, but behave as permanently methylated or demethylated when occupied, respectively, by aspartate (D) or glutamine (Q). In addition, substitution ofh þ 3 by asparagine (N) in mutant E(N)EEE is also informative as it partially impairs the methylation of site 3, which would correspond to a partial mehthylation state in our model.  (0) Normalized methylation rates for different CheB − mutants reproduced from refs. 28,30 . The four letters (D for aspartate, E for glutamate, N for aspargine, and Q for glutamine) in the first column are the residues, respectively, in the methylation sites 1, 2, 3, and 4. The four middle columns are the methylation rates of each site in arbitrary unit. Simulations mimicking each mutant were performed with their configurations shown in the last column, with x representing modifiable sites. We fixedμ 4 ¼ 0 due to the low methylation rate of site 4.  We studied the methylation dynamics of these CheB − mutants by using a sequential-dominant model with a small value of η (= 0.1). As shown in Fig. 6a, the sequential methylation of a completely demethylated receptor (0(1)xxx(0)) begins by methylating site 3, afterwards proceeds to site 2, and to site 1. This order, i.e., p 3 (t) > p 2 (t) > p 1 (t), persists throughout the methylation process consistent with experiments results. We also studied the methylation dynamics when the starting siteμ 3 is fixed to beμ 3 ¼ 0 andμ 3 ¼ 1 to mimic the EEDE and EEQE receptors, respectively. As shown in Fig. 6b, when we fixμ 3 ¼ 1, the order of methylation for site 2 and site 1 still persists, i.e., p 2 (t) > p 1 (t), which is again consistent with experiments.
Finally, the most informative and also most stringent test of our theory comes from the mutant receptor EEDE. Besides a much slower methylation rate, an inverted behavior p 2 < p 1 was observed experimentally in EEDE (Table 1). Remarkably, these behaviors in particular the inversion also appear in our model as shown in Fig. 6c. The reason for this inversion is that sequential modification is broken when site 3, the starting site in the sequence, cannot be methylated. As a result, the downstream sites (site 2 and site 1) have to be methylated (at least initially) by the random methylation process, which has no a priori preference between site 1 and site 2. Once site 2 becomes methylated, it will enhance the methylation rate at site 1 due to the sequential methylation process but not the other way around. Thus the sequentiality between site 2 and site 1 leads to the observed inversion. Consistent with this argument and with the role played by the initiator, the partial methylation ofh þ 3 in EEE(N)E reduces the ratio p 3 /p 2 , when compared to EEEE. Quantitatively, the CheB − data lead to an estimate for a lower bound for the random methylation parameter η ≥ 0.047 (see Supplementary Discussion for more details), which is consistent with the estimated range of η for Tar from the wt data above.
Overall, our model results, together with existing experimental data, suggest that the methylation process for sites 3, 2, and 1 are mostly sequential and are affected by the initiatorh þ 3 , but there is a small but finite random component.
Testable predictions for future experiments. Our model can be used to predict the methylation level profile for different mutants, which can be tested by future experiments. As the reference, the methylation levels of the wild-type cell (0(1)xxx(0) receptors in the presence of CheR and CheB) decrease monotonically from site 3 to 1 as shown in Fig. 7a. We first study the mutant withh þ 3 ¼ 0, which inhibits sequential methylation of site 3. As shown in Fig. 7b,h þ 3 ¼ 0 brings down the methylation of site 3, leading to site 2 being more methylated than sites 1 and 3. To explore the inhibition of sequential demethylation of site 1, we next study the mutant withh À 1 ¼ 1. As shown in Fig. 7c, site 2 is less methylated than sites 1 and 3 in the steady state. Finally, we study the mutant withh À 1 ¼ 1 andh þ 3 ¼ 0 in which both methylation of site 3 and demethylation of site 1 are inhibited. As shown in Fig. 7d, the steady state methylation profile monotonically increases from 3 to 1, which is exactly the inverse of the wt profile (Fig. 7a).
We can also predict the effects of mutating the key methylation sites (μ 1 andμ 3 ) on adaptation dynamics. We first studied effects of mutating site 1 or site 3 to be permanently unmethylated by fixing either μ 3 = 0 or μ 1 = 0 in our model. We found that adaptation still works in μ 3 = 0 mutant [0(1)xx0(0)] but is severely impaired in the μ 1 = 0 mutant [0(1)0xx(0)] as shown in the Supplementary Fig. 7a. We next studied effects of mutating site 1 or site 3 to be permanently methylated by fixing either μ 3 = 1 or μ 1 = 1 in our model. We found that response to decrease in attractant concentration remains intact in theμ 3 ¼ 1 mutant [1(1)1xx(0)], but is severely impaired in theμ 1 ¼ 1 mutant [1(1)xx1(0)] as shown in Supplementary Fig. 7b. These predictions can be tested by measuring the kinase activity dynamics in vivo in these mutants by using FRET 32 .

Discussion
Multisite regulatory proteins are ubiquitous in biology, yet their functions are not well understood. Here, we studied effects of ordering among the multiple modification sites and possible benefits of having multiple sites in the context of bacterial chemotaxis. We discuss the two main findings below.
First, we found that sequential modification is crucial for perfect adaptation. Previous study 14,35 showed that perfect adaptation can be achieved by an integral control mechanism where dynamics of the controller (receptor methylation level) only depend on receptor activity. Here, we showed that sequential modification is another important ingredient for the integral control mechanism as it guarantees the methylation/demethylation rates to be independent of the receptor methylation level, Eq. (4). As a direct consequence of sequential modification, the adapted activity is independent of the receptor methylation level (or the stimulus strength), i.e., perfect adaptation.
We note that there may be other possible scenario for the methylation/demethylation rates to be independent of the available modification sites. For bacterial chemoreceptors, the binding and unbinding of CheR to the receptor are faster than its catalytic rate and the dissociation constant K D is relatively small 36 . If the enzyme binds to all available active sites randomly with equal probability, the number of available sites effectively changes the substrate concentration. Given that the substrate concentration is much higher than the Michaelis-Menten constant K M ≈ K D , the methylation reaction rate, which is limited by the slow catalytic reaction, would be independent of the substrate concentration and thus independent of the number of available methylation (demethylation) sites. However, the binding rate (k on ) of the enzyme in this scenario would depend on the substrate concentration and the available modification sites, which seems to be inconsistent with the recent in vitro measurements of the k on rates for CheR binding to Tar(EEEE) and Tar(QQQQ) receptors 36 (see Supplementary Note 1 for more details). Furthermore, the random methylation pattern predicted by this scenario is inconsistent with the observed sequentiality among the different methylation sites in in vivo experiments 28,30 .
Second, we found that there is a tradeoff between response gain and adaptation accuracy. We showed that this tradeoff can be improved significantly by having more modification sites but only with the sequential modification process. Taken together, our study suggests a general two-pronged strategy to enhance chemotaxis performance by having multiple modification sites to extend the dynamic range of high gain, and a sequential modification process to maintain adaptation accuracy. Direct comparison with existing experiments confirms our theory and reveals that the methylation process for methylation sites 3, 2, and 1 of Tar is mostly sequential with a small but finite random component 0.05 ≤ η ≤ 0.13. The confirmed importance of sequential receptor methylation begs the question of the underlying molecular mechanism responsible for maintaining specific ordering in multisite modification. Previous mutant studies showed that methylation of a given site is affected by a residue seven amino acids to the C terminus 28,30 , which is exactly how sites 1, 2, and 3 are arranged (Fig. 5). Also, methylation of site 3 is affected by a residue 7 amino acids to the C terminus, even though that residue itself is not a methylation site 28 . Indirect evidence of sequential demethylation by cheB can also be found in refs. 37,38 (see Supplementary Discussion for details).
The existing experiments mentioned above suggest a chain reaction scheme for the sequential methylation process. However, it is not clear whether the preceding site in the sequence increases the binding affinity of CheR to the receptor or the catalytic rate or both. It is also not clear whether and how different receptors in the closely packed receptor cluster compete for the limited CheR molecules in the cluster. We believe that a detailed biochemical model that incorporates key steps such as binding/unbinding and catalytic reactions in the methylation/demethylation processes together with quantitative in vitro measurements of the methylation/demethylation rates for wt and mutant receptors are needed to address these questions. The same strategy should also be used to study the much less known demethylation process. In addition to searching for possible molecular mechanisms for ordered modification, another interesting question is what are the thermodynamic costs of implementing such ordered modification mechanisms for accurate control 24 . Finally, it is worth pointing out that even though the detailed molecular mechanism of the methylation and demethylation reactions still remains open, our conclusions regarding the general properties of the system such as response gain, adaptation accuracy, and their tradeoff and their dependence on the level of sequentiality (η) of the underlying multisite modification process should hold true.
Our work serves as a successful case study of multisite protein modification by using a modeling approach in combination with knowledge of the underlying biochemical pathway and quantitative data. This combined approach provides a powerful general framework that can be applied to other signaling systems to understand the mechanisms of multisite signaling proteins and their biological functions.

Methods
The standard model for bacterial chemotaxis. We briefly describe a previously developed general mathematical framework-the standard model for studying bacterial chemotaxis signaling pathway dynamics (see ref. 17 for a recent review).
In the standard model for bacterial chemotaxis, each receptor has two key state variables-its kinase activity (a) and its methylation state (μ). For kinase activity, a receptor can be either active (a = 1) or inactive (a = 0). For methylation state, as each receptor has M(≥1) modification (methylation) sites, there are a total of 2 M possible modification states characterized by a M − dimensional binary vector μ = (μ 1 , μ 2 , … , μ M ), where the binary number μ j = 0, 1 respectively represents the unmethylated and methylated state of site j(=1, 2, … , M). The total modification level of a receptor is given by: m k μ k ¼ P M j¼1 μ j . The receptor kinase activity dynamics is fast relative to its methylation dynamics. Here, we use the standard two-state model to describe the receptor kinase activity dynamics, where the active and inactive states are separated by a free energy difference Δf. When the fast ligand-receptor binding/unbinding process is averaged out, Δf(m, [L]) depends on the receptor's total modification level m and the ligand concentration [L]. From previous studies on bacterial chemotaxis 18,35 , the free energy Δf can be written as: where K I and K A are the dissociation constants of the ligand binding to the inactive and active conformations of the receptor, α(<0) is the free energy change due to adding (or removing) one methylation group to the receptor, and m 0 determines the average modification level in the absence of any stimulus ([L] = 0). Another important phenomenon in bacterial chemotaxis is that bacterial chemoreceptors form polar clusters [39][40][41] . The receptors and their kinase activities are coupled with each other in the cluster. Following previous work 19,42 , we model the receptor cluster by using an Ising-type model with nearest neighbor interaction with strength C.
Dynamic Monte-Carlo simulations of the Ising-type model (see Supplementary Methods for details of the Monte-Carlo simulations) is used to obtain the distribution of receptors, P aμ , in a given state (aμ), which describes the statistical properties of the receptor cluster. From the full distribution function P aμ , distribution of the microscopic methylation state μ can be obtained by summing over the fast variable a, P μ ¼ P 1 a¼0 P aμ , and the probability of the total modification level m is given as: From these distribution functions, average properties of the receptor cluster can be obtained. For example, the average methylation level is According to Eq. (10), kinase activity of a receptor a h i m only depends on its total methylation level m, which can be expressed as: and the average activity for all receptors is: These distribution functions and average receptor properties are used here to understand the response gain and adaptation accuracy in bacterial chemotaxis quantitatively. In particular, we focus on investigating how different modification schemes (random or sequential) affect the adaptation accuracy and response gain in this paper.
Characterizing the performance of the chemotaxis signaling pathway. The performance of the chemotaxis signaling pathway can be characterized by two key system-level properties: the integrated response gain (amplification) Γ and adaptation accuracy ξ −1 , which we define in the following.
At a given background ligand concentration [L], the adapted methylation levels of all receptors in the system (receptor cluster) are represented by m A ð½LÞ (vector m = (m 1 , m 2 , … , m N ) contains the methylation levels of all the receptors in the receptor cluster, m i is the methylation level of receptor-i (1 ≤ i ≤ N) and N is the number of receptors in the cluster) and the average adapted activity of the system is given by: a h i A ð½LÞ a h iðm A ð½LÞ; ½LÞ. Upon a sudden change of ligand concentration from [L] to [L] + δ[L], the system first responds by a change of activity δ〈a〉, which can be written as: where the negative sign is due to the fact that increase of attractant concentration leads to decrease of receptor activity in bacterial chemotaxis. To describe the system's ability to amplify the input stimulus over a broad range of background stimulus concentration, we define the overall gain Γ as the integral To characterize the adaptation accuracy over the wide range of backgrounds, we define an overall adaptation error ξ by integrating ϵ over the stimulus concentration in log-scale (natural base is used here for convenience): The overall adaptation accuracy is defined as the inverse of the adaptation error, ξ −1 . Details of the Monte-Carlo simulations and of computing Γ and ξ are discussed in Supplementary Methods.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All data used to support the findings of this work are available upon request.

Code availability
The code used to perform the simulations is available at https://github.com/ bernardomello/chemotaxis.