Main

Large 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).

Spontaneous two-state switching in chemosensory arrays

To 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.

Switching temporal statistics collapse to those of two-dimensional Ising model

To 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…N**j, where N**j 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 cooperativity

We 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 < 0.6kBT), on rescaling by factors (\epsilon =\left|;{J}^{* }-J\right|/J) and L−(z−b) according to the finite-size scaling theory (Extended Data Fig. 7 shows the determination of scaling parameters). Right margin: estimated distribution (Supplementary Fig. 8) of the switching timescale ratio (\left\langle \Delta t\right\rangle /\left\langle \tau \right\rangle) across cells, for Tsr-I214K (red) or Tar [QEEE] (blue). The shaded areas (left margin) correspond to one standard deviation of the estimated timescale ratio distribution. b, Phase diagram of the switching behaviour as in the L–J plane. Parameter value pairs generated polarized (grey region) and non-polarized (dark blue region) fluctuations as estimated from the peak-valley ratio of the time series histograms. Isolines of similar switching timescale ratios, as indicated, were constructed from the simulations (beaded coloured and black curves). Theoretical finite-size scaling predictions for the critical coupling energy J*(L) (ref. 62; gold curve) and the switching time ratios (dashed lines; Supplementary Equation (9) and Supplementary Table 16 list the fit parameters) are superimposed. Inset: switching timescale ratio isolines for Tar (blue) and Tsr (red), normalized by the finite-size scaling of the critical coupling energy J* reveals values within 3% of the critical energy.

Because different combinations of J and L can yield the same value of ϵ**L, this scaling does not allow the unique determination of J from the measured values of r (Fig. 3a). However, finite-size scaling theory predicts that the critical coupling energy J*, separating the highly ordered (polarized) and disordered (non-polarized) regimes, also depends on the system size L (refs. 61,62). Consistent with this prediction, a phase diagram on the J–L plane constructed from simulations showed that the boundary between polarized (Fig. 3b, grey region) and non-polarized (Fig. 3b, blue region) dynamics coincided well with the theoretically predicted scaling (Fig. 3b, gold curve) J*(L) ≈ J*(∞)/(1 − c**L−1), where J*(∞) ≈ 0.44kBT is the critical coupling energy in the thermodynamic limit (L → ∞) and c is a boundary-condition-dependent constant (c = 1.25 for free boundary conditions62; Supplementary Note 2.3). Interestingly, we found that isolines ({J}_{r}^{{,{\rm{iso}}}(L)) corresponding to L–J combinations that yield the same timescale ratio r (Fig. 3b, beaded curves) were nearly parallel with the profile of J*(L). Dividing J*(L) from the isolines ({J}_{r}}{,{\rm{iso}}}(L)) corresponding to the measured r values for Tar (rTar ≈ 9; Fig. 3b, blue) and Tsr (rTsr ≈ 12; Fig. 3b, red) revealed a remarkable confinement of the coupling energy for both Tar and Tsr to within ±3% of J*(L) across a broad range in L (Fig. 3b, inset). Thus, despite uncertainty in the array size L, finite-size scaling analysis of the observed temporal statistics strongly suggests that both Tar and Tsr arrays are poised very close to the Ising critical point.

Arrays without adaptation feedback exhibit critical slowing of response

Finite-size scaling analysis also allowed us to estimate the fundamental flipping timescale 1/ω0 of allosteric units. Combining the observed residence and transition times with an approximate range for the size L based on the known structure, stoichiometry and component expression levels of the array, the scaling could be calibrated to yield 1/ω0 ≈ 10–35 ms (Extended Data Fig. 8 and Supplementary Note 3). This allowed us to directly compare the response dynamics between simulations and experiments. Previous theoretical studies have shown that signal amplification due to coupling within Ising-type conformational spread models comes at the cost of prolonged response times47,63, potentially impairing temporal comparisons during run-and-tumble navigation64,65,66. To address the influence of near-critical cooperativity on response times to external stimuli, we first simulated the dynamic response of arrays with J ≈ J* after an applied step stimulus (implemented within the model by an external ligand field ΔH that modifies the flipping rate ω0 by a factor of eΔH(2a−1)/2; Methods) to mimic kinase inactivation on attractant chemoeffector stimulation. These simulations demonstrated that near-critical Ising arrays indeed respond slowly, with a broad distribution of times averaging tens of seconds (Fig. 4a). We then conducted FRET experiments with a sub-saturating stimulation of adaptation-deficient cells expressing homogeneous arrays (consisting of a single chemoreceptor species, Tsr) using a polydimethylsiloxane (PDMS) microfluidic flow cell that allowed for chemoeffector concentration changes much faster than the acquisition frequency of 1 Hz (Methods). The experimental r