Spontaneous switching in a protein signalling array reveals near-critical cooperativity

MainLarge protein assemblies at the heart of many cell signalling processes exhibit varying degrees of structural and dynamical order, from liquid-like granules1,2,3 to solid-like arrays4,5,6. Recent experiments are revealing how phase transitions that lead to the formation of liquid-like assemblies—discontinuous phase transitions analogous to the condensation of gas into liquid—are used by cells to implement various information processing tasks, such as signal initiation and confinement7, kinetic proofreading8 and noise control9. Theory also predicts that different phase transitions can occur within solid-like assemblies—continuous phase transitions, analogous to the ordering of magnetic spins in a ferromagnet at low temperatures—arising from conformational interactions between protein subunits10. However, requirements for continuous phase transitions are more stringent than those for discontinuous phase transitions—they occur only at a special point in phase space—a critical point. It remains unclear if such critical transitions occur in signalling assemblies and, if so, how they constrain or contribute to the functional design of protein assemblies. Elucidating the physical principles of signalling assemblies would not only advance our understanding of cell signalling in nature but also provide new design principles for the rational design of synthetic protein circuits11,12,13,14.Recent structural studies are revealing an increasing number of solid-like signalling assemblies with a high degree of spatial order, where subunit proteins are arrayed in a regular pattern. The repertoire of such crystal-like assemblies reported so far is diverse in both form and function, and includes one-dimensional filament-like assemblies found in cellular homeostasis and inflammation signalling15,16, protein rings that mediate the control of cell motility and apoptosis17,18, as well as two-dimensional arrays involved in chemosensing, neuromuscular control and innate immune responses19,20,21,22. Yet, despite exciting advances in resolving the ultrastructure of these large assemblies23,24, the mechanistic design principles of their signalling function remain challenging to address experimentally. Of particular interest is the ability of signalling assemblies to perform signal processing by way of cooperative interactions between assembly subunits. In contrast to the compact oligomeric signalling complexes of fixed size, these extensive structures tend to assemble through open-ended polymerization and/or multivalent interactions, making them inherently variable in both size and composition5. The resulting polydispersity and stoichiometric diversity tend to mask their true signalling dynamics25,26, as well as their size and composition dependencies6,27, both in vivo and in vitro. An ideal functional assay would, therefore, interrogate assembly level dynamics in singulo—at the level of an individual assembly—but experimental realization has remained elusive.Here we report in vivo FRET experiments that resolve this challenge for the chemosensory array of Escherichia coli, a canonical two-dimensional signalling assembly lining the cytoplasmic membrane, which allows motile bacteria to bias their random-walk swimming trajectories to navigate chemical environments28,29,30. This higher-order assembly comprises thousands of receptor, kinase and scaffolding molecules arranged in a well-defined lattice structure31,32,33 and, thus, serves as a paradigm for signalling in two-dimensional protein assemblies. The array integrates input signals (chemoeffector ligand concentrations) with a negative feedback signal (covalent modification state of receptors) to generate a signal output (activity of the kinase, CheA) that can be followed in real time by an intermolecular Förster resonance energy transfer (FRET) system28. By labelling two downstream proteins, the response regulator CheY, which is phosphorylated by CheA, and its phosphatase CheZ, which dephosphorylates CheY, a FRET read-out proportional to the kinase activity can be obtained within live cells28 (Fig. 1a). Our strategy leverages recent advances in extending this in vivo FRET technique to the single-cell level34,35, to achieve functional measurements of individual arrays within live cells. Although this FRET system provides a whole-cell measurement that integrates the kinase output of all arrays within each cell, the number of arrays per cell is highly variable36 and in a substantial minority of cells, nearly all receptors and other array-component proteins are assembled within one dominant large array37. Thus, by searching for cells in which signalling is dominated by a single large array, we aimed to achieve in singulo FRET measurements of functional array output.Fig. 1: In singulo measurements of chemosensory array dynamics reveal two-state switching fluctuations.a, Illustration of single-cell FRET assay to measure the kinase activity34, with a schematic of the side view of a membrane-associated chemosensory array of chemoreceptors (either Tar or Tsr). The kinase (CheA) within the chemosensory array phosphorylates CheY (CheY-P), whereas CheZ dephosphorylates CheY. The FRET signal (black arrows) between fluorescently labelled CheZ and CheY, measured from the emitted fluorescence intensities (red and yellow arrows) is proportional to the activity of the chemosensory array. b,c, Activity time series from single-cell FRET experiments that reveal two-state switching in the absence of any chemoeffector, of representative cells without adaptation enzymes expressing only the chemoreceptor Tar [QEEE] (b; blue) or the chemoreceptor Tsr-I214K (c; red). Cells are exposed to the measurement buffer (MotM) for most of the experiment. To obtain the minimum and maximum FRET levels, cells are exposed to an attractant stimulus (Tar: 1 mM α-methylaspartate (MeAsp); Tsr: 1 mM l-serine (Ser), grey bands) and a repellent stimulus (Tar: 0.3 mM NiCl2; Tsr: 1 mM l-leucine (Leu), purple bands). These levels were then used to normalize the activity time series of each cell between zero and unity. Example time series other than the two-state one are shown in Extended Data Fig. 1. d, Chemosensory cluster organization from CheZ localization (left) and the corresponding activity time series (right, raw in grey, 6 s (averaged); in colour) and histograms (on the right margin) of activity determined by FRET, for a representative cell expressing only the chemoreceptor Tsr-I214K, with a single visible cluster exhibiting two-state switching (orange) and a representative cell with two visible clusters exhibiting multistate switching (green). Additional time series and cluster organization are shown in Extended Data Fig. 2. e, Histograms of the number of states in switching cells, for cells with one visible cluster (top, orange) and for cells with two visible clusters (bottom, green).Full size imageSpontaneous two-state switching in chemosensory arraysTo identify cells whose chemoreceptor population is concentrated into a single large array, we focused on fluctuations in the kinase output that are observable by FRET at the single-cell level. In wild-type (WT) Escherichia coli cells, the stochastic kinetics of the reversible covalent modification reactions (methylation/demethylation) mediated by the adaptation system are known to drive temporal fluctuations in kinase activity under constant environmental conditions34,35,38. Surprisingly, however, previous single-cell FRET studies identified the largest fluctuations in engineered genetic mutants devoid of the adaptation system and expressing only one of the five chemoreceptor species34,35, with a subpopulation of cells exhibiting two-level (all-ON/all-OFF) fluctuations in kinase output as the input ligand signal was held constant, at a sub-saturating level34 (Supplementary Note 1). If these fluctuations were driven by the intrinsic dynamics of the array, such synchronized two-level stochastic switching of the entire kinase population would be expected only if whole-cell signalling were dominated by a single dominant array. However, because an external ligand was present in those experiments, the possibility remained that the observed output switches were caused by spurious fluctuations in the external chemoeffector signal. We therefore began by testing whether two-level fluctuations occur in the absence of any external ligand stimulus by exploiting known genetic modifications of chemoreceptors that mimic the effects of methylation/demethylation, and can modify the chemoreceptor activity bias and maintain otherwise normal signalling function39. We tuned down the activity bias by expressing the aspartate receptor Tar in the QEEE covalent modification state (from QEQE in WT Tar), which yields an intermediate activity bias without the addition of chemoeffectors (Supplementary Note 1 and Extended Data Fig. 1e)39. This allowed us to measure temporal fluctuations by single-cell FRET recordings in thousands of individual cells (Fig. 1b and Extended Data Fig. 1a). In 204 out of 1,414 cells obtained from 19 FRET experiments, we observed two-level switching between a high- and low-activity state (>65% of transitions showing activity-level changes of >0.7). Spontaneous switching in the absence of exogenous ligands was not specific to Tar receptors, as we observed a similar switching behaviour (two-state switching in 548 out of 4,446 cells across 44 experiments) in analogous experiments with Tsr-I214K, a single-residue replacement mutant of the serine receptor Tsr within its ‘control cable’ region40 with a downshifted activity bias similar to that of Tar [QEEE] (Fig. 1c and Extended Data Fig. 1b).Spontaneous two-state switching is a hallmark of allosteric signalling complexes such as ion channels41, but their observation requires measurements at the level of individual complexes; ensemble-averaged experiments cannot resolve the switching dynamics because the timings of switching events are uncorrelated and their dynamics are averaged out. To further test whether the observed two-level switches represent the behaviour of individual chemosensory arrays, we performed additional single-cell FRET experiments at a reduced expression level of the CheY/CheZ FRET pair. The attenuated cytoplasmic fluorescence in these experiments allowed the detection of chemosensory arrays as the intensity peaks of CheZ-YFP fluorescence, due to the well-established phenomenon of CheZ clustering at chemosensory arrays42,43, and the number of stable kinase output states could be determined through FRET experiments on the same individual cell (Fig. 1d and Extended Data Fig. 2). We found that cells that exhibited only a single detectable cluster predominantly demonstrated only one or two stable activity states, whereas cells with two clusters typically demonstrated three or more states (Fig. 1d and Extended Data Fig. 2). Tests with mutants deficient in either ligand binding or response cooperativity, as well as tests under a metabolic perturbation, further supported the idea that the observed switches are driven by intrinsic stochastic dynamics of the array, rather than extrinsic fluctuations (Extended Data Figs. 3 and 4 and Supplementary Note 1). Collectively, these results strongly support the idea that two-level fluctuations observed in our FRET experiments reflect cooperative signalling dynamics within a single dominant array.To consider the physical mechanism underlying this long-range cooperativity, we analysed the temporal statistics of switching fluctuations (Fig. 2a) by extracting the time interval between switching events, hereafter called the residence times Δtup and Δtdown for time spent in the up (a ≈ 1) and down (a ≈ 0) states, respectively, as well as the duration of the activity transient on switching, hereafter called the transition times τ+ and τ− for upward and downward switches, respectively (Fig. 2a). We first interpreted these data as a barrier-crossing stochastic process in which the up and down states correspond to wells within an energy landscape (Fig. 2b), the shape of which could be approximated from the observed activity time series histograms with the free energy difference ΔG (in units of the thermal energy kBT) between the up and down states determined by the activity bias \(\left\langle a\right\rangle\) as \(\Delta G={\mathrm{ln}}[(1-\left\langle a\right\rangle )/\!\left\langle a\right\rangle ]\). Consistent with this, we found that for both Tar [QEEE] and Tsr-I214K arrays, the residence time intervals between switching events were exponentially distributed across the full range of \(\left\langle a\right\rangle\) (Fig. 2c), and the average residence time of each cell \(\left\langle \Delta {t}_{{\rm{up,down}}}\right\rangle\) as a function of ΔG was found to obey an Arrhenius-type exponential scaling \(\left\langle \Delta {t}_{\mathrm{up,down}}\right\rangle (\Delta G)=\langle \Delta t\rangle {{\rm{e}}}^{-{\gamma }_{\mathrm{up,down}}\Delta G}\) (Fig. 2d), where 〈Δt〉 is a characteristic residence timescale independent of the cell’s activity bias and γup,down are fitting constants. From the crossings of the Arrhenius-fit lines in Fig. 2d, we determined 〈Δt〉 for both receptor species: 〈Δt〉Tar = 47.0 ± 1 s and 〈Δt〉Tsr = 65.5 ± 1 s.Fig. 2: Temporal statistics of switching events are well described by a two-dimensional Ising model.a, Definitions of residence times Δtup,down and transition times τ+,−, determined routinely for each switching event (Methods). b, Coarse-grained energy landscape along the array-activity coordinate a (estimated as the negative logarithm of the activity histogram) based on the selected time series with 〈a〉 ≈ 0.5 (15 cells; left) and 〈a〉 ≈ 0.9 (12 cells; right) from cells expressing Tsr-I214K. Horizontal scale bar indicates the reaction coordinate Δa = 1, vertical scale bar shows the energy (in kBT) and the dashed lines indicate the free energy difference ΔG between the high- and low-activity state. c, Histogram of residence times from experiments. Residence times for Tsr-I214K (left) and Tar [QEEE] (right), with each event sorted by the activity bias of the corresponding cell (colours as in d). In each panel, data (points) are shown together with fits to single exponential functions (solid lines). Fit parameters and number of data points are shown in Supplementary Tables 3, 5, 12 and 13. pdf, probability density function. d, Mean residence times per cell as a function of the energy bias ΔG between the high- and low-activity state for all cells expressing Tsr-I214K (top) and Tar [QEEE] (bottom). Fit parameters and number of data points are shown in Supplementary Table 7. e, Conformational spread model of the chemosensory array activity with size L × L, in which each individual unit can switch between the activity states of 1 (white) or 0 (dark). A difference in neighbouring spins is associated with an energy cost of J, shown for three different transitions on a lattice with size L = 4 in the absence of an external field (H = 0). f, Example activity time series obtained by simulating dynamics on a strongly coupled (L = 26, J = 0.4625kBT, dark grey) and weakly coupled (L = 20, J = 0.2375kBT, light grey) lattice. The strongly coupled lattice exhibits stochastic switches between two activity levels (dashed lines). The simulated time series was downsampled to approximate the experimental acquisition frequency (solid black line). g, Same data as in c, but histograms of the residence times from the simulated time series with L = 12 and J = 0.5kBT, sorted by the activity bias generated by an applied external field H. Fit parameters and number of data points are shown in Supplementary Tables 8, 10, 14 and 15. h, Histograms of transition times for Tsr-I214K (solid lines) and Tar [QEEE] (dashed lines). Mean transition times and number of data points are shown in Supplementary Tables 4 and 6. i, Transition times from simulated two-state time series (N = 12, J = 0.5kBT, varying H). To approximate the experimental signal-to-noise ratio (Supplementary Fig. 6), Gaussian white noise was added to the simulated time series. Mean simulated transition times and number of data points are shown in Supplementary Tables 9 and 11. Experimental (points) and simulated (solid lines) histograms are scaled horizontally to have mean transition time of cells expressing Tsr, for comparison. Error bars represent 95% confidence intervals obtained through bootstrap resampling.Full size imageSwitching temporal statistics collapse to those of two-dimensional Ising modelTo investigate whether and how the observed temporal statistics could be explained by cooperative subunit interactions, we turned to theory. We used an Ising-type conformational spread model of allosteric cooperativity44,45,46,47, which assumes that subunit conformations are coupled through nearest-neighbour interactions. The strength of these interactions is parameterized by a coupling energy J (promoting order) whose magnitude relative to kBT (promoting disorder) determines a finite spatial range (that is, a correlation length) over which action at one site can affect distant sites. By varying J, therefore, the model can represent allosteric systems along a continuous scale of conformational disorder, including that of the more widely used Monod–Wyman–Changeux model (which is recovered on taking the fully ordered limit J → ∞ and has been applied extensively in modelling chemosensory arrays48,49,50,51). Importantly, for signalling function, two-dimensional Ising models are known to exhibit a continuous phase transition as a function of J, with a spontaneous ordering of subunit conformations above a critical coupling energy J*. Although so far, strong experimental support for conformational spread models has been obtained in one-dimensional protein rings52,53, Ising models exhibit a critical point only in two or higher dimensions, and the implications of the Ising phase transition in two-dimensional protein assemblies remained untested experimentally.We performed kinetic Monte Carlo simulations on an L × L lattice of allosteric units with free boundary conditions, each of which can flip between two conformational states, active (a = 1) or inactive (a = 0) (Fig. 2e). The flipping rate of the unit at site i was modified from a fundamental flipping frequency ω0 by the influence of its nearest neighbours (j ∈ 1…Nj, where Nj is the number of nearest neighbours) through the coupling energy J (in units of kBT) as \(\omega ={\omega }_{0}\exp[{-J(2{a}_{i}-1){\sum }_{j}(2a_{j}-1)}]\), corresponding to an Ising model in which each active–inactive bond on the lattice contributes an energy penalty of J (Methods). At a low coupling strength (J ≪ J*), each unit switches independently and the total activity of the array demonstrates only small fluctuations about its mean value at \(\left\langle a\right\rangle =1/2\) (Fig. 2e). However, as the coupling energy is increased towards its critical value J*, the correlation length approaches the finite size of the array, generating a double-well potential and the arrays exhibit switching events between fully active and fully inactive states (Fig. 2f and Extended Data Fig. 5). We analysed the temporal activity statistics of a simulated array with parameters within this two-state switching regime, with various values of a weak biasing field Hb that modifies the flipping rate by a factor \({{\rm{e}}}^{{H}_{{\rm{b}}}(a-1/2)}\), to approximate the diverse FRET activity biases observed across individual cells in the population (Extended Data Fig. 1e,f). Simulated residence time distributions (Fig. 2g and Extended Data Fig. 5d) were in excellent agreement with their experimental counterparts (Fig. 2c), recapitulating their characteristic exponential shape at each activity bias.By contrast, the measured transition time distributions had peaked profiles for both Tar and Tsr arrays (Fig. 2h). The average downward transition time \(\left\langle {\tau }_{-}\right\rangle\) and the average upward transition time \(\left\langle {\tau }_{+}\right\rangle\) were similar between Tar and Tsr arrays, with \(\left\langle {\tau }_{-}\right\rangle\) slightly greater than \(\left\langle {\tau }_{+}\right\rangle\) in both cases (\(\left\langle {\tau }_{+}^{\,{\rm{Tsr}}}\right\rangle =4.29\pm 0.06\) s, \(\left\langle {\tau }_{-}^{\,{\rm{Tsr}}}\right\rangle =6.07\pm 0.07\) s, \(\left\langle {\tau }_{+}^{\,{\rm{Tar}}}\right\rangle =4.79\pm 0.08\) s, and \(\left\langle {\tau }_{-}^{\,{\rm{Tar}}}\right\rangle =6.06\pm 0.09\) s; mean ± s.e.m.). These modest yet significant differences between \(\left\langle {\tau }_{+}\right\rangle\) and \(\left\langle {\tau }_{-}\right\rangle\) hint at the breaking of time-reversal symmetry and possible non-equilibrium driving54,55, and are not captured by the equilibrium Ising model. Remarkably, however, when normalized by their respective mean values to remove this asymmetry, we observed a remarkable collapse of all measured transition time distributions onto the profile of the simulated distributions (Fig. 2i). Furthermore, both measured and simulated transition-time statistics demonstrated no dependency on the activity bias (Extended Data Figs. 5e,f and 6). Collectively, the high degree of quantitative agreement between these measured and simulated temporal statistics suggest that an Ising-type conformational spread model with near-critical coupling strength (J ≈ J*) provides an excellent approximation to chemoreceptor array dynamics.Finite-size scaling analysis reveals near-critical cooperativityWe sought to quantify the degree to which both Tar and Tsr arrays are close to criticality. Crucial in considering critical phenomena in living systems are finite-size effects56, which modify the quantitative behaviour near critical points compared with well-known results derived in the thermodynamic limit (where the system size L → ∞). Finite-size scaling theory of Ising-type models is highly developed57,58, but direct comparisons between the temporal statistics of the Ising model and our experimental data are complicated by the fact that the fundamental flipping timescale 1/ω0 of cooperative units (corresponding to the spin-flip timescale in the Ising model) remains unknown. We, therefore, identified—as a key experimental observable—the dimensionless ratio \(r\equiv \left\langle \Delta t\right\rangle /\left\langle \tau \right\rangle\) between the residence and transition timescales (with \(\left\langle \tau \right\rangle\) defined as the average \(\left\langle \tau \right\rangle \equiv (\left\langle {\tau }_{+}\right\rangle +\left\langle {\tau }_{-}\right\rangle )/2\); Supplementary Note 2), which divides out the ω0-dependence to enable direct comparisons. Simulations at various values of J indeed revealed a strong dependence of this ratio r on L (Fig. 3a), with all results collapsing onto a single curve defined by the finite-size scaling relation \(r\approx {L}^{(z-b)}\exp ({c}_{0}\epsilon L)\), where z, b and c0 are scaling constants and \(\epsilon =\left|\;{J}^{* }-J\right |/J\) is the ‘reduced temperature’ providing a dimensionless measure of the (energetic) distance to criticality59,60 (Supplementary Note 2.5 and Extended Data Fig. 7 detail the determination of the scaling constants).Fig. 3: Finite-size scaling analysis reveals near-critical cooperativity of chemosensory arrays.a, Left: switching timescale ratio \(r=\left\langle \Delta t\right\rangle/\left\langle \tau \right\rangle\) for various values of the coupling energy J (blue to black) as a function of lattice size L (circles), with exponential fits (dashed lines). Inset: data collapse for the near-critical region (J 10 s), limiting the speed at which bacteria can respond. Taken together, these results provide a clear experimental demonstration that cooperative interactions within the array affect the signalling response timescale of E. coli. Although cooperativity-induced response slowdown is most acute in cells engineered to express only a single chemoreceptor species, it is substantial also in cells with the WT complement of chemoreceptor species.Near-critical cooperativity balances a speed–amplitude trade-offWhy are bacterial chemosensory arrays poised so close to criticality, despite potentially deleterious slowing of response? To address this question, we examined the relationship between response speed and response amplitude using simulations for various combinations of J and ΔH (Fig. 5a,b). For every stimulus size ΔH, increasing the coupling strength J led to a decrease in the response speed (defined as the inverse of response time) but an increase in response amplitude, indicating a trade-off (Fig. 5c). The profile of these HL isolines are interesting when viewed as a Pareto front70 for navigating the trade-off, having a convex shape with a knee above which the response speed drops off sharply. Remarkably, the critical coupling energy (J = J*; Fig. 5c, gold curve) traverses this knee at every stimulus size, indicating that near-critical cooperativity enables a balancing of these two response objectives, allowing for large response amplitudes without drastically compromising the response speed.Fig. 5: Near-critical cooperativity balances response amplitude and speed.a, Top: without adaptation feedback, near-critical cooperativity extends the range of cooperativity across the entire array, leading to all-or-none response. Shown are example states (binary squares) of an array initialized in the all-active (\(\left\langle a\right\rangle \approx 1\)) state at time t = 0 when a positive external field (ΔH > 0, mimicking chemoattractant addition) was applied, and at a later time t > tR exceeding the response time tR of the array. b, Simulated response time series (dashed curves) of activity bias \(\left\langle a\right\rangle (t)\), computed by averaging 48 stochastic trajectories of a 20 × 20 array without adaptation feedback and J = 0.47kBT, for various values of added external field ΔH (c shows the colour code), initialized at 〈a〉 ≈ 1 and assuming 1/ω0 = 30 ms (Extended Data Fig. 8). Inset: example single-array stochastic trajectories (solid curves), which invariably demonstrate all-or-none switching-like response. The coloured rectangles illustrate the definition of the response time tR as the average time of crossing half-maximal activity, a = 0.5 (coloured rectangle). The response amplitude is defined as \(1-\left\langle a\right\rangle (J,\Delta H)/0.5\). c, Speed–amplitude trade-off in array responses without adaptation feedback. For each stimulus size ΔH (see the legend for colour code), a thin curve indicates the dependence of response amplitude and speed on the coupling energy J. The thick gold curve is the critical isoline, tracing out points on each thin curve corresponding to \(J={J}_{J}^{* }\). The horizontal dashed line represents the speed required to respond within a typical run time (1 s) of E. coli92, corresponding to an approximate lower bound for effective run-and-tumble chemotaxis. d, Adaptation feedback injects spatial disorder within the array that limits the extent of cooperativity and stabilizes intermediate-activity states. Shown are example array states (binary squares) in the presence of adaptation feedback (indicated by grey arrows), before stimulus (t tA), each with corresponding activity bias \(\left\langle a\right\rangle\). e, Simulated response time series (dashed curves) of activity bias 〈a〉(t), computed by averaging 96 stochastic trajectories of a 20 × 20 array with adaptation feedback and J = 0.47kBT, for various values of added external field ΔH (f shows the colour code), initialized at 〈a〉 ≈ 0.5. Response time and amplitude are determined from an exponential fit to the averaged activity time series (coloured rectangles). Inset: example single-array stochastic trajectories (solid curves), smoothed with a 7-s moving-average filter to aid visual inspection. f, Same data as c, but for speed–amplitude trade-off in array responses with adaptation feedback. Note that response slowdown is mitigated compared with non-adapting arrays, and the speed of near-critical arrays (thick gold curve) remains above the lower-bound response speed set by the run-and-tumble behaviour timescale (horizontal dashed line) across all the stimulus sizes.Full size imageHowever, these Ising simulations also indicate that at J = J*, the speed of response to small-amplitude stimuli (H ={\int }_{-\infty }^{\infty }{\int }_{-l/2}^{+l/2}I(x,y)\,{\rm{d}}x\,{\rm{d}}y,$$ (3) Then, the cluster size is defined as$$\Delta I={I}_{{\rm{m}}{\rm{a}}{\rm{x}}}- ,$$ (4) where Imax is the maximum of the fluorescence intensity along the long axis of the cell I.Switching analysisSwitching events were analysed automatically using a custom-made MATLAB script (Mathworks, version 2019b or newer). For each cell, attractant and repellent responses were detected automatically based on the timing and duration of the stimulus delivery. Next, the ratiometric FRET time series of each cell, after correcting for photobleaching, was normalized between one and zero based on the repellent and attractant response amplitudes, respectively. The FRET signal was then low-pass filtered using a 3-s moving-average filter, and switching events were detected as peaks in the derivative of the filtered signal. With each switching event, an amplitude (change in kinase activity), residence time (time until next event) and transition time (duration of the switch based on a fit of the form 1 − e−1) are associated. The switching behaviour of each cell was classified according to the amplitude and number of the switching events per cell. A two-state switching cell is defined as a cell with at least 65% of its transitions showing activity-level changes of at least 0.7 or 70% of total kinase activity—increasing these thresholds to 80% and 0.8, respectively, has only a marginal effect on the switching statistics (Supplementary Fig. 5). For each two-state switching cell, the bleaching correction is refined by another correction step using only the activity in either the a = 0 state, for cells with bias 0.5. The time series is then renormalized using the maximum activity level of the histograms during the two-state switching, and the amplitude, transition time and residence time associated with each switching event are extracted again. To verify that the automated switching analysis extracts events reliably, we generated mock time series that recapitulated the switching phenotype of Tsr-I214K and added Gaussian white noise to approximate the experimental signal-to-noise ratio. By comparing the transition and residence times of the mock time series to the times extracted by our automated analysis, we find that the relative uncertainty is minimal, even for low signal-to-noise ratios (Supplementary Fig. 6). Events occurred at a stable average frequency throughout the course of experiments (Supplementary Fig. 6). There were no strong correlations between the transition and residence times per cell (Supplementary Fig. 7). The variability between cells was estimated by convolving simulated transition and residence time distributions with Gaussian white noise to match the experimentally observed distributions (Supplementary Fig. 8).Numerical simulations of the conformational spread modelConformational spread was modelled by a two-dimensional Ising model on an L × L lattice with free boundary conditions. Each lattice site represents an allosteric unit, whose conformational state was represented in the main text by an activity variable ai ∈ {0, 1}. For the discussion here, we make the mapping σ = 2a − 1 and discuss the model in terms of the ‘spin variable’ σ, which takes one of two values: σ = 1 for active, and σ = −1 for inactive. The activity of the unit at site i is influenced by its Nj-nearest neighbours through the coupling energy J a biasing field Hb and a ligand field HL, giving a Hamiltonian (total energy) \({\mathcal{H}}\) for this lattice (in units of kBT):$${\mathcal{H}}=-J\mathop{\sum }\limits_{\left\langle ij\right\rangle }{\sigma }_{i}{\sigma }_{j}+({H}_{{\rm{b}}}/2+{H}_{{\rm{L}}}/2)\mathop{\sum }\limits_{i}{\sigma }_{i},$$ (5) where \(\left\langle ij\right\rangle\) indicates that the summation is over all the nearest-neighbour pairs in the lattice. The biasing and ligand fields (Hb and HL, respectively) were set to zero for most of our simulations, except where specified otherwise in the text (the field dependence of the switching phenotype is described in more detail in Supplementary Note 2.7). The probability of finding the lattice in a given configuration is then proportional to \({{\rm{e}}}^{-{\mathscr{H}}}\) and the ratio of probabilities pi(σi) and pi(−σi) for the ith site to be in state σi as opposed to −σi is$$\frac{{p}_{i}({\sigma }_{i})}{{p}_{i}(-{\sigma }_{i})}={{\rm{e}}}^{-\Delta {\mathcal{H}}},$$ (6) where \(\Delta {\mathcal{H}}\equiv {\mathcal{H}}({\sigma }_{i})-{\mathcal{H}}(-{\sigma }_{i})\) is the change in the Hamiltonian on flipping σi. We set the rate ω(σi → −σi) for the unit at site i to flip from a state σi to the opposite state −σi as$$\omega ({\sigma }_{i}\to -{\sigma }_{i})={\omega }_{0}\exp \left[-J{\sigma }_{i}\mathop{\sum }\limits_{j}^{{N}_{j}}{\sigma }_{j}+({H}_{{\rm{b}}}/2+{H}_{{\rm{L}}}/2){\sigma }_{i}\right]$$ (7) where ω0 is the fundamental flipping frequency of a single lattice unit and Nj is the number of nearest neighbours, so as to satisfy (together with equations (5) and (6)) the detailed balancing condition:$$\frac{{p}_{i}({\sigma }_{i})}{{p}_{i}(-{\sigma }_{i})}=\frac{\omega (-{\sigma }_{i}\to {\sigma }_{i})}{\omega ({\sigma }_{i}\to -{\sigma }_{i})}.$$ (8) We note that the fundamental frequency ω0 corresponds to the rate of conformational transitions in an individual allosteric unit. Because this rate is unknown for the bacterial chemosensory array, in our simulations, we set this parameter to unity, implying that the simulated temporal statistics are expressed in units of the fundamental timescale 1/ω0.The rates, as defined in equation (7), were used in a kinetic Monte Carlo scheme (essentially as in ref. 101, but modified to have free boundary conditions, as shown below), which uses one random number in every iteration to draw the time until the next flip, and another to determine which site flips. The lattice was an L × L lattice with free boundary conditions, such that the number of nearest neighbours Nj at each lattice site was$${N}_{j}=\left\{\begin{array}{ll}2 & \,{\rm{at\; corners,}}\\ 3 & \,{\rm{at\; edges,}}\\ 4 & \,{\rm{otherwise.}}\end{array}\right.$$ (9) To calculate the activity a of the lattice at each time point, we map the spin variable σi ∈ {−1, +1} at each site to an activity variable ai ∈ {0, 1} as ai = (σi + 1)/2 and take its mean across the lattice:$$a=\frac{1}{{L}^{2}}\mathop{\sum }\limits_{i}{a}_{i}.$$ (10) To simulate adaptation feedback, a bias field is defined in which each lattice unit has a specific bias field Hb,i that represents the effect of receptor methylation according to39$${H}_{{\rm{b}},i}=\alpha ({m}_{i}-{m}_{0}),$$ (11) where 0 0) are methylated with rate kR. This activity-dependent (de)methylation activity results in robust perfect adaptation to a steady-state activity level of a0 = kR/(kB + kR). We defined these rates in terms of the fundamental frequency ω0 according to kB = kR = nR/L2 × VR/ω0 = 0.0015ω0, where nR = 200 is the approximate number of adaptation proteins per cell103 and VR = 1/10s is the maximum rate of the methylation reaction104, reflecting that the adaptation enzymes work at saturation105. Arrays were simulated for a sufficient time to reach the steady state (a0 ≈ 0.5) before the application of an external field HL (representing chemoattractant ligand stimulation) was used to probe the array responses.Simulations were implemented with Python 3.7 or newer, and single simulation runs for extended times were performed on regular desktop computers. Parallel computations were performed on the LISA cluster of the SURFsara national computing facility (Amsterdam). For the parallel runs, each parameter set was given a unique seed for random number generation. The resulting time series were sampled at regular intervals, and subsequently downsampled to approximate the acquisition frequency of experiments (1 Hz) relative to the array-level switching frequency (~10−2 Hz). Simulated activity time series were then further processed in the same way as experimental data to extract the temporal statistics of switching.Calibration of fundamental frequencyAlong the isoline \({J}_{r}^{\,{\rm{iso}}}(L)\) (Fig. 3b), the timescale ratio r by definition remains constant, meaning 〈Δt〉 = r〈τ〉. Combined with the calibration of the finite-size scaling for the transition time 〈τ〉sim ≈ Lb through numerical simulations (Extended Data Fig. 7), and noting that timescales from experiments and simulations are related as \({\langle \Delta t\rangle }_{\exp }={\omega }_{0}^{-1}{\langle \Delta t\rangle }_{{\rm{sim}}}\), we obtain ω0 = rcτLb/Δtexp, where cτ and b are constants based on the scaling analysis (Supplementary Note 2). The resulting dependence on L was very similar for both Tar and Tsr (Extended Data Fig. 8). Although our FRET measurements do not provide a direct estimate of L, we can motivate the approximate upper and lower bounds based on structural and biochemical findings in the literature between L = 17 and L = 30 (Supplementary Note 3). With these approximate limits, our scaling analysis yields a flipping timescale of individual allosteric units in the range 1/ω0 ≈ 15–35 ms (Extended Data Fig. 8).Determination of correlation lengthThe correlation length \({\mathcal{E}}(J)\) was determined from exponential fits to the correlation function of Ising arrays106:$$c(r)=\overline{{\sigma }_{i}{\sigma }_{j}}={{\rm{e}}}^{\frac{-r}{{\mathcal{E}}(\,J)}},$$ (12) where \(\overline{{\sigma }_{i}{\sigma }_{j}}\) is the spin product averaged over all pairs with equal distance r = ∣i − j∣. For each unit i, j, pairs are sampled horizontally and vertically (for example, keeping either i or j constant, ignoring diagonal values). For each J, the correlation length \({\mathcal{E}}(J)\) is computed as the average correlation length of 48 independent array states.Statistics and reproducibilityNumber of experiments and number of cells per experiments are described in the figure captions or in supplementary tables. Two-state switching cells were selected according to the criteria described above; otherwise, no data were excluded from the analysis. No method was used to predetermine the sample sizes. The investigators were not blinded to allocation during experiments and outcome assessment.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.