# Signature of (anti)cooperativity in the stochastic fluctuations of small systems: application to the bacterial flagellar motor

John Palmeri, Maria-José Franco-Oñate, Nils-Ole Walliser, Andrea Parmeggiani, Jérôme Dorignac, Frédéric Geniet, Jean-Charles Walter, Francesco Pedaci, Ashley Nord

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

John Palmeri, Maria-José Franco-Oñate, Nils-Ole Walliser, Andrea Parmeggiani, Jérôme Dorignac, et al.. Signature of (anti)cooperativity in the stochastic fluctuations of small systems: application to the bacterial flagellar motor. International Physics of Living Systems (IPOLS) 2023, Daniel Goldman, Georgia Tech, Aug 2023, Atlanta, United States. [hal-04890444](https://hal.umontpellier.fr/hal-04890444v1)

HAL Id: hal-04890444

https://hal.umontpellier.fr/hal-04890444v1

Submitted on 29 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Signature of (anti)cooperativity in the stochastic fluctuations of small systems: application to the bacterial flagellar motor

María-José Franco-Oñate¹, Andrea Parmeggiani¹, Jérôme Dorignac¹, Frédéric Geniet¹, Jean-Charles Walter¹, Francesco Pedaci², Ashley L Nord², John Palmeri¹,*, and Nils-Ole Walliser¹,*

¹Laboratoire Charles Coulomb (L2C), Univ. Montpellier, CNRS, Montpellier, France

²Centre de Biologie Structurale (CBS), Univ. Montpellier, CNRS, INSERM, Montpellier, France

*Corresponding authors: john.palmeri@umontpellier.fr, nils-ole.walliser@umontpellier.fr

## Abstract

The cooperative binding of molecular agents onto a substrate is pervasive in living systems. To study whether a system shows cooperativity, one can rely on a fluctuation analysis of quantities such as the number of substrate-bound units and the residence time in an occupancy state. Since the relative standard deviation from the statistical mean monotonically decreases with the number of binding sites, these techniques are only suitable for small enough systems, such as those implicated in stochastic processes inside cells. Here, we present a general-purpose grand canonical Hamiltonian description of a small one-dimensional (1D) lattice gas with either nearest-neighbor or long-range interactions as prototypical examples of cooperativity-influenced adsorption processes. First, we elucidate how the strength and sign of the interaction potential between neighboring bound particles on the lattice determine the intensity of the fluctuations of the mean occupancy. We then employ this relationship to compare the theoretical predictions of our model to data from single molecule experiments on bacterial flagellar motors (BFM) of *E. coli*. In this way, we find evidence that cooperativity controls the mechano-sensitive dynamical assembly of the torque-generating units, the so-called stator units, onto the BFM. Furthermore, in an attempt to quantify fluctuations and the *adaptability* of the BFM, we estimate the stator-stator interaction potential. Finally, we conclude that the system resides in a *sweet spot* of the parameter space (*phase diagram*) suitable for a smoothly adaptive system while minimizing fluctuations.

## Contents

1 Introduction 2

2 Results 3

2.1 Lattice gas model 3

2.2 Mean occupancy and its standard deviation 3

2.2.1 System behavior at half-filling 5

2.3 Fluctuations of the mean occupancy and their dependency on $\langle\varphi\rangle$ and $J$ 6

2.3.1 Fluctuations in the zero-cooperativity limit: the Hill-Langmuir case 7

2.3.2 Fluctuations in the strong cooperativity limit 8

2.3.3 Fluctuations for intermediate finite cooperativity 8

2.3.4 Fluctuations in the strong anti-cooperativity limit 9

2.4 Probability distribution function (PDF) of the occupancy and its dependency on $\langle\varphi\rangle$ and $J$ 9

2.4.1 PDF in the zero-cooperativity limit: the Hill-Langmuir case 10

2.4.2 PDF in the strong cooperativity limit 10
<table><tr><td>2.4.3</td><td>PDF in the strong anti-cooperativity limit</td><td>11</td></tr><tr><td>2.4.4</td><td>PDF near half-filling: summary</td><td>11</td></tr><tr><td>2.5</td><td>Application to the bacterial flagellar motor (BFM)</td><td>11</td></tr><tr><td><b>3</b></td><td><b>Discussion and Conclusion</b></td><td><b>14</b></td></tr><tr><td><b>4</b></td><td><b>Methods</b></td><td><b>16</b></td></tr><tr><td>4.1</td><td>1D short-range lattice gas (SRLG)</td><td>16</td></tr><tr><td>4.2</td><td>Probability distribution function of the occupancy</td><td>18</td></tr><tr><td>4.3</td><td>Infinite range (IR) model</td><td>23</td></tr><tr><td>4.4</td><td>System size, fluctuations, and experimental precision</td><td>25</td></tr><tr><td>4.5</td><td>Effective Hill coefficient</td><td>26</td></tr><tr><td>4.6</td><td>Supplemental plots</td><td>27</td></tr><tr><td>4.7</td><td>Experimental data</td><td>30</td></tr></table>

# 1 Introduction

Cooperativity is a pervasive phenomenon that emerges in biological systems at different time and length scales, and levels of complexity [1]: from the molecular scale [2–5] to multiple-cell organization [6], extending to the emergence of collective behavior in many-body systems.

We focus here on cooperative processes involving the adsorption of ligands onto a substrate disposing of a limited number of binding sites [7]. Our main goal is to develop a general method for assessing how short- and long-range interactions between substrate-bound ligands affect stochastic fluctuations in the number of adsorbed units. In the presence of cooperativity, small system size effects strongly influence the amplitude of the stochastic fluctuations and the shape of the equilibrium probability distribution function (PDF) describing substrate occupancy. By studying these characteristic signatures of cooperativity, we propose 1) a criterion to determine whether any given adsorption system exhibits cooperative or anti-cooperative behavior and 2) a method to quantify the amplitude of the ligand-ligand interaction potential.

To this end, we introduce a minimal 1D lattice gas model with interacting units in equilibrium with a thermal bath and chemical reservoir, described using a general-purpose grand canonical Hamiltonian with short or long-range interactions. With it, we study how fluctuations at equilibrium, described by the standard deviation of the occupancy, are influenced by the model parameters, namely the ligand binding energy, the ligand-ligand interaction potential, and the chemical potential of the bulk reservoir. To get closer to observable quantities, we invert the problem by determining the key formulae relating the occupancy standard deviation and PDF to the experimentally accessible average occupancy of the system. In this way, we are able to quantify how, at fixed average occupancy, increasing cooperativity leads to an increase in fluctuations and how increasing anti-cooperativity leads to a decrease. In the process, we highlight the saturation of fluctuations for both strong cooperativity and strong anti-cooperativity, as well as a subtle difference between fluctuations near half-filling for systems with even and odd binding sites in the case of strong anti-cooperativity. Finally, we apply the model by comparing our theoretical predictions with experimental data for the bacterial flagellar motor (BFM) of *E. coli* [8, 9], a macro-molecular complex that has previously been found to exhibit cooperative binding of the torque-generating (stator) units [10]. In this previous work, although evidence for cooperativity was uncovered in stator unit assembly dynamics, no estimate of the stator-stator interaction potential was proposed, which is our goal here. In Methods Section 4.5, we make the connection between our model-based approach and phenomenological approaches based on effective Hill coefficients [7, 11].

We also discuss the concept of motor adaptability to changes in external conditions and set up a general framework to assess its desired level depending on how the motor is expected to
function. Our fluctuation analysis provides evidence for cooperativity by leading to an estimate of the stator-stator interaction potential between 1 and 2 $k_BT$ for the short-range model, a result that is coherent with this system's expected smooth adaptive nature to external stimuli [12–16]. We contrast this moderate cooperative behavior found for stator unit binding to the BFM with the highly cooperative behavior of rotor-switching units that are believed to be behind the rapid stochastic switching in BFM rotational direction [11, 17].

# 2 Results

In this section, we propose a minimal prototypical model of cooperative particle binding on a substrate of finite size. We then obtain an analytical expression for the fluctuations of the mean occupancy as a function of the relevant model parameters, including the interaction potential between bound particles, which acts as a proxy for system cooperativity. Finally, we apply these results to investigate the characteristic signatures of cooperativity by establishing a general criterion for its manifestation and then propose a procedure to estimate the interaction potential from experimental data.

## 2.1 Lattice gas model

Consider a periodic 1D lattice with $L$ binding sites in contact with a heat and particle reservoir of constant temperature $T$ and effective chemical potential $\mu$. We assign to each binding site $i$ ($i = 1, ..., L$) a binary variable $\phi_i$: $\phi_i = 1$ means that a particle (ligand) occupies the $i$-th site; otherwise, $\phi_i = 0$. The array $\phi = \{\phi_i\} = (\phi_1, ..., \phi_L)$ uniquely defines one of the $2^L$ possible microscopic configurations (or microstates). The relative occupancy in a given microstate is $\phi = L^{-1} \sum_{i=1}^L \phi_i$. We say that the system is in mesostate $N$ if the number of particles on the lattice, i. e. the occupancy of the system equals $N = L\phi$. There are $C_N^L = L!/[N!(L-N)!]$ microstates $\phi$ consistent with the occupancy $N$; in this case, we say that the mesostate $N$ has multiplicity $C_N^L$.

The energy of the system depends on the specific microscopic configuration $\phi$ according to the grand canonical effective 1D short-range interaction Hamiltonian:

$$ \beta\mathcal{H}(\phi) = -J \sum_{i=1}^{L} \phi_i \phi_{i+1} - \mu \sum_{i=1}^{L} \phi_i. \qquad (1) $$

Since we choose periodic boundary conditions, $i = L + 1$ corresponds to $i = 1$. Furthermore, we introduce $\beta = 1/(k_B T)$ with the Boltzmann constant $k_B$, rendering the two adjustable control parameters of our model, $J$ and $\mu$, dimensionless. The effective chemical potential, $\mu$, is considered to be made up of two parts, $\mu = \mu_r - \epsilon$, where $\mu_r$ is the chemical potential of the external particle reservoir and $\epsilon < 0$ is the binding energy of particle adsorption onto the substrate. The nearest-neighbor interaction $J$ is attractive for $J > 0$ (cooperativity) and repulsive for negative values (anti-cooperativity). For fixed $J$ the number of bound particles can be modulated by varying $\mu$ through $\mu_r$ and/or $\epsilon$. When $J = 0$ we recover the usual non-interacting (Hill-Langmuir) model.

The predictions for this short-range model will also be compared below with those for the long-range model where all pairs of particles interact with the same strength independent of their relative distance on the lattice.

## 2.2 Mean occupancy and its standard deviation

The Hamiltonian Eq. (1) falls into the universality class of the 1D short-range lattice gas model. Hence, we can map our model to the 1D Ising model (see [18] for a comprehensive overview). It is, therefore, possible to derive analytical expressions for the grand partition function and
<figure><img src="image_3.png" /><figcaption>Figure 1: Contour plots of the mean relative occupancy in equilibrium $\langle\phi\rangle$ (left), and the standard deviation of the relative occupancy in equilibrium $\sigma = \sqrt{\langle\phi^2\rangle - \langle\phi\rangle^2}$ (right), as functions of the dimensionless interaction potential $J$ and chemical potential $\mu$, according to expressions from Eq. (2) and (4), for a system of size $L = 13$.</figcaption></figure>

both the first and second moments of the relative occupancy in thermodynamic equilibrium by employing the transfer matrix formalism. By taking the derivative of the grand partition function $\Xi$ with respect to $\mu$, we can calculate the mean relative occupancy at equilibrium, $\langle\phi\rangle = L^{-1} \sum_{i=1}^{L} \langle\phi_i\rangle = (L\Xi)^{-1} \partial_{\mu}\Xi$:

$$ \langle\phi\rangle = \frac{1}{2} \left[ 1 + \tanh \left( \frac{L}{2\xi} \right) \frac{\sinh X}{\sqrt{\sinh^2 X + e^{-J}}} \right], \quad (2) $$

where $X = (J + \mu)/2$, and $\xi$ is the equilibrium correlation length of the system (given in number of lattice sites), which can be rewritten in terms of $J$ and $\mu$ as follows (we refer to Section 4.1 for more details):

$$ \xi = 1 / \ln \left( \frac{\cosh X + \sqrt{\sinh^2 X + e^{-J}}}{\cosh X - \sqrt{\sinh^2 X + e^{-J}}} \right). \quad (3) $$

The mean equilibrium occupancy equals $\langle N \rangle = L\langle\phi\rangle$. For purposes of illustration in what follows, unless otherwise stated, we shall use $L = 13$. This particular choice will allow us to compare with greater ease the theoretical results presented here to the experimental data when we later, in Section 2.5, apply our model to the study of the recruitment of torque-generating units (or stator units) onto the BFM. It should be kept in mind that our theoretical results are valid for any value of $L$ and that our conclusions for small systems apply to any sufficiently small value of $L$.

We propose that the operational definition of “small system” be simply those systems for which fluctuations, as measured by the standard deviation of relative occupancy, are within experimental resolution and can be used to extract the sign and amplitude of the interaction strength within the usual biophysical range. A more careful rendering of this definition can be found in Section 4.4. Such a definition is not contingent on the amplitude of (anti)cooperativity and therefore can encompass even the non-interacting (Hill-Langmuir) model.

Fig. 1a depicts the mean relative occupancy from Eq. (2) as a function of the two control parameters of the system: $J$ and $\mu$. The half-filling (HF) contour line $\langle\phi\rangle = 0.5$ cuts the plot diagonally, meaning that the lattice is, on average, half-filled whenever $J = -\mu$ (or $X = 0$) holds. Below (above) that line, the system’s mean occupancy is lower (higher) than 50 %.
<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison between the average relative occupancy and the standard deviation as functions of $\mu$ for two values of cooperativity: $J = 0$ (Hill-Langmuir) and $J = 2$.</figcaption></figure>

Contour lines run together close to the plot's bottom right sector diagonal, corresponding to positive $J$ and negative $\mu$. In this sector defined by $X \approx 0$, for $J \gg 1$, the mean occupancy for a fixed positive value of the interaction potential becomes extremely sensitive to slight chemical potential variations. Systems characterized by parameters in that range behave like biological switches, nearly jumping from zero to full occupancy with small changes in $\mu$. On the contrary, more adaptive systems, that is, systems whose occupancy varies more smoothly from zero to full occupancy as a function of $\mu$, would be located in the central and upper left sector of Fig. 1a.

For the short-range lattice gas model from Eq. (1), the standard deviation $\sigma$ of the occupancy at equilibrium can be expressed in terms of the derivative of the mean relative occupancy with respect to the chemical potential as follows (cf. Section 4.1): $\sigma = \sqrt{\langle\phi^2\rangle - \langle\phi\rangle^2} = L^{-1/2}\sqrt{\partial_\mu\langle\phi\rangle}$. Applying this result to Eq. (2), we obtain a compact analytical expression for $\sigma$:

$$ \sigma = \frac{1}{2} \left[ \text{sech}^2\left(\frac{L}{2\xi}\right) \frac{\sinh^2 X}{\sinh^2 X + e^{-J}} + \tanh\left(\frac{L}{2\xi}\right) \frac{e^{-J} \cosh X}{L(\sinh^2 X + e^{-J})^{3/2}} \right]^{1/2}, \quad (4) $$

whose dependency on $J$ and $\mu$ is depicted in Fig. 1b. In what follows, we will use $\sigma$ as a measure of fluctuations, hence the stochastic nature of the system.

In Fig. 2, we compare the average relative occupancy, $\langle\phi\rangle$, and the standard deviation, $\sigma$, as functions of $\mu$ for two vertical cuts in the above contour plots: $J = 0$ (Hill-Langmuir) and $J = 2$. We note that with increasing $J$ the $\langle\phi\rangle$ curve sharpens and shifts to lower values of $\mu$ and that the maximum in $\sigma$ increases in amplitude and also shifts to lower values of $\mu$. These shifts are due to the displacement with increasing $J$ of the value of $\mu$ at half-filling, $\mu_{HF}(J) = -J$, where the slope of the $\langle\phi\rangle$ curve and the amplitude of the standard deviation are both largest.

### 2.2.1 System behavior at half-filling

Fluctuations increase as $J > 0$ grows. Along the half-filling diagonal ($X = 0$) in Fig. 1b, they saturate at a finite global maximum as $J$ goes to infinity.

By evaluating Eq. (4) and (3) for $X = 0$, one obtains an explicit formula for the standard deviation at half-filling (HF) as a function of $L$ and $J$:

$$ \sigma_{\text{HF}}(J) = \frac{e^{J/4}}{2\sqrt{L}} \sqrt{\tanh\left(\frac{L}{2\xi_{\text{HF}}(J)}\right)}, \quad (5) $$

with

$$ \xi_{\text{HF}} = 1/\ln\left(\frac{1+e^{-J/2}}{1-e^{-J/2}}\right) = 1/\ln[\coth(J/4)], \quad (6) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Parametric plots of the standard deviation of the relative occupancy, $\sigma$, versus the mean relative occupancy, $\langle\varphi\rangle$, for different values of the interaction potential: from $J = -10$ (bottom curve) to $J = 10$ (upper curve) in steps of 1. Plots (a) and (b) correspond to systems of size $L = 13$ and $L = 14$, respectively.</figcaption></figure>

which represents the maximum equilibrium correlation length at fixed $J$ for any value of mean occupancy.

In the Hill-Langmuir limit ($J = 0$), the correlation length vanishes, and the standard deviation at HF becomes $1/(2\sqrt{L})$. For low couplings ($0 \le J < 2$), the correlation length increases sub-linearly with $J$ and is accurately given by $\xi_{HF} \approx -1/\ln(J/4)$. In this same limit, the standard deviation increases exponentially with $J$, $\sigma_{HF}(J)/\sigma_{HF}(0) \approx e^{J/4}$, an approximation valid for $J < 4$, when $L \approx 10$. This increase can be interpreted as a widening of the probability distribution function of states of a given $N$ around the average value $\langle N \rangle$, where the probability is a maximum. In the high cooperativity limit $J \gg 1$, $\xi_{HF}$ shows an asymptotic exponential growth, $\xi_{HF} \sim \frac{1}{2}e^{J/2}$, which is the analog of the well-known zero temperature critical behavior of spin chains [19]; and, consequently, for fixed $L$, $\sigma_{HF}$ tends to $1/2$ when the correlation length surpasses $L/2$ and the system evolves into a bimodal one (for which only the zero and full occupancy states contribute, see below). We remark that although for fixed $J$, $\sigma_{HF}$ vanishes as $L^{-1/2}$ when $L \to \infty$, as expected, in the limit where $L$ is kept fixed and $J$ tends to $+\infty$, the standard deviation at HF reaches the finite (maximum, or saturation) value of one-half.

Fluctuations decrease as $J < 0$ becomes more negative. Along the half-filling diagonal ($X = 0$) in Fig. 1b ($L = 13$), they saturate at a finite value as $J \to -\infty$. This is a general result for systems with an odd number of sites. In this strong anti-cooperative limit, for small enough systems, it turns out that fluctuations depend sensitively on whether the number of lattice sites is even or odd. For $L$ even, $\sigma_{HF}$ tends to zero; for odd $L$, it tends to $1/(2L)$. This leads to the striking result that in the limit of strong anti-cooperativity, a fluctuation analysis near HF could clearly detect the difference between small systems with even and odd numbers of sites, provided that the experimental precision is high enough.

## 2.3 Fluctuations of the mean occupancy and their dependency on $\langle\varphi\rangle$ and $J$

Often, the experimental parameter one can measure (or control indirectly, as is the case for the BFM by varying the external load, cf. Section 2.5) is not the chemical potential of the reservoir $\mu$, but the mean number of bound ligands at equilibrium, $\langle N \rangle$. It is, therefore, insightful to plot the standard deviation $\sigma$ as a function of the mean relative occupancy $\langle\varphi\rangle$, which, in general, can be done using a parametric plot.

To unfold this relation, we express $\langle\varphi\rangle$ and $\sigma$, for different but fixed values of $J$, as functions of $\mu$ and plot this family of curves in Fig. 3a and Fig. 3b for an odd and even number of lattice sites, respectively. Each curve for a given $J$ represents the set of points of coordinates $(\langle\varphi\rangle = f_J(\mu), \sigma = g_J(\mu))$ with $\mu$ spanning from -15 to 30, which turns out to be a suitable one-dimensional domain to complete the plots. We have chosen the interaction potential $J$ to take values from -10 to 10 in steps of size 1. The standard deviation shows a left-right
particle-hole symmetry with respect to the system half-filling: $\langle\varphi\rangle_{HF} = f_J(-J) = 1/2$, which is a manifestation of the underlining particle-hole symmetry of the system under study. The system has zero fluctuations when $\langle\varphi\rangle = 0$ (empty substrate) or $\langle\varphi\rangle = 1$ (full substrate), which is easy to understand because a distribution of states with non-zero probability around these average values would be incompatible with the average values themselves. For non-negative values of $J$, the maxima of fluctuations appear at half-filling, taking on the value of $1/(2\sqrt{L})$ for $J = 0$ and $1/2$ for $J \to +\infty$, and the $\sigma$ vs. $\langle\varphi\rangle$ curves converge asymptotically towards the limiting $J \to +\infty$ one. More subtly, results for the anti-cooperative limit $J \to -\infty$ depend on whether $L$ is odd or even.

A key result of the above analysis and one of our major conclusions is that for cooperative systems of size $L \approx 10$ and mean relative occupancies not too close to 0 or 1, in the range $0 < J < 7$ the curves of $\sigma$ vs. $\langle\varphi\rangle$ are spaced sufficiently far enough apart that the standard deviation is suited for estimating $J$. The same conclusion can be drawn for anti-cooperative systems of size $L \approx 10$ and mean relative occupancies close to $1/2$ in the range $-3 < J < 0$. These results are illustrated in Fig. 3a for $L = 13$ and Fig. 3b for $L = 14$. Outside these ranges, the standard deviation saturates to limiting values, and therefore experimental data near these limits could only be used to put bounds on the value of $J$.

We will investigate the cooperative and anti-cooperative regimes in more detail in what follows, taking care to study those limits for which we were able to derive simplified exact or approximate analytical expressions for $\sigma(\langle\varphi\rangle)$.

### 2.3.1 Fluctuations in the zero-cooperativity limit: the Hill-Langmuir case

In the absence of cooperativity ($J = 0$), which corresponds to the case of a simple adsorption process with on-site volume exclusion interactions only, the free energy calculated from Eq. (1) reduces to an expression determined entirely by the mesostate (the occupancy) $N$: $\beta H_0(\varphi|N) = -\mu N$. The correlation length from Eq. (3) vanishes; hence each adsorbing site acts as an independent particle trap. Thus we recover the Hill-Langmuir adsorption model at equilibrium for which $\langle\varphi\rangle$ takes the following typical sigmoid form as a function of $\mu$:

$$ \langle\varphi\rangle_0(\mu) = 1/(1+e^{-\mu}). \qquad (7) $$

The corresponding standard deviation becomes:

$$ \sigma_0(\mu) = \frac{\operatorname{sech}(\mu/2)}{2\sqrt{L}} = \sqrt{\frac{\langle\varphi\rangle_0 - \langle\varphi\rangle_0^2}{L}}, \qquad (8) $$

which can be reinterpreted to obtain

$$ \sigma_0(\langle\varphi\rangle) = \sqrt{\frac{\langle\varphi\rangle - \langle\varphi\rangle^2}{L}}. \qquad (9) $$

This latter expression exemplifies the general result that, in the thermodynamic limit, $L \to \infty$, $\langle\varphi^2\rangle \to \langle\varphi\rangle^2$, and hence, as expected, the standard deviation (of the relative occupancy) tends to zero as $L^{-1/2}$ (self-averaging property at equilibrium in the thermodynamic limit).

In Section 4.2, we show using a Gaussian approximation for the PDF that for $J = 0$ the standard deviation is inversely proportional to the curvature (in absolute value) of the entropy of mixing, which reaches a minimum at HF. Therefore, the standard deviation for $J = 0$ reaches a maximum there. More generally, one can show that for $J > 0$ not too large, a Gaussian approximation for the PDF is accurate and therefore the standard deviation is inversely proportional to the curvature (in absolute value) of an effective free energy, which reaches a minimum at HF.
### 2.3.2 Fluctuations in the strong cooperativity limit

Away from HF, although $\xi$ remains finite in the large $J$ limit, it can still be large compared with $L/2$ for not-too-large values of $L$, provided that one is not too close to zero or full filling (see Fig. 15). Taking the limits $J \to \infty$ and $\mu = 2X - J \to -\infty$ simultaneously, but keeping $L$ finite and $X$ constant, allows Eqs. (2) to be simplified to

$$ \langle\varphi\rangle_{\infty}(X) = \frac{1}{2} \left[ 1 + \tanh \left( \frac{L}{2\xi_{\infty}(X)} \right) \operatorname{sgn}(X) \right], \quad (10) $$

and $\langle\varphi^2\rangle_{\infty} = \langle\varphi\rangle_{\infty}$, where $\xi_{\infty}(X) = 1/(2|X|)$, leading to an infinite $J$ standard deviation of

$$ \sigma_{\infty}(X) = \frac{1}{2} \operatorname{sech} \left( \frac{L}{2\xi_{\infty}(X)} \right) = \sqrt{\langle\varphi\rangle_{\infty} - \langle\varphi\rangle_{\infty}^2}, \quad (11) $$

which can be reinterpreted to obtain

$$ \sigma_{\infty}(\langle\varphi\rangle) = \sqrt{\langle\varphi\rangle - \langle\varphi\rangle^2}. \quad (12) $$

The maximal standard deviation at half-filling is due to the large width of the distribution of states with non-zero probability compatible with $\langle\varphi\rangle = 1/2$ (see Section 2.4). This limit is easy to understand in terms of a bimodal system: only the empty and full states contribute because the probability of intermediate states is suppressed by the presence of energetically costly domain walls (as will be discussed in detail below). The bimodal system, therefore, behaves as an effective Hill-Langmuir model with a system size equal to 1, as can be seen by taking $L=1$ in Eqs. (9).

### 2.3.3 Fluctuations for intermediate finite cooperativity

We have seen in Fig. 3 how, at fixed average occupancy, increasing cooperativity leads to an increase in fluctuations. A simple, insightful, and accurate way of understanding this increase can be formulated by starting from the previously derived result for the Hill-Langmuir model and introducing the concept of block domains. When $J=0$, all sites are decoupled, and the standard deviation is non-zero for a finite-size system simply because of statistical fluctuations (absence of self-averaging). When $J>0$, neighboring sites are coupled and can be grouped (approximately) into correlated block domains of size $b(\xi, L) > 1$, leading to

$$ \sigma(J; \xi, L) \approx L^{-1/2} \sqrt{\langle\varphi\rangle - \langle\varphi\rangle^2}, \quad (13) $$

where $L_{\text{eff}} \equiv L/b$. When $1 < 2\xi < L$, correlations lead to block domains that we take to be of size $b=2\xi$. The above approximation allows the system to be described by a reduced number, $L_{\text{eff}}$, of independently fluctuating block domains of size $b$. Since $b$ should also tend to 1, when $\xi \to 0$ (zero cooperativity, Hill-Langmuir); and to $L$, when $2\xi/L \gg 1$ (strong cooperativity, bimodal behavior), a convenient interpolation formula between the above three cases is

$$ b(\xi, L) = \frac{\tanh[L/(2\xi)]}{\tanh[1/(2\xi)]}. \quad (14) $$

By comparing the predictions of this approximation with the exact results at half-filling, where the standard deviation is a maximum, one can see that this approximation is extremely accurate (see Fig. 14).

Exact results for the standard deviation can be obtained by writing the mean square relative occupancy at equilibrium,

$$ \langle\varphi^2\rangle = \frac{1}{L^2\Xi} \partial_{\mu}^2 \Xi = L^{-2} \sum_{i=1}^{L} \sum_{j=1}^{L} \langle\varphi_i\varphi_j\rangle, \quad (15) $$
in terms of the 2-point correlation function, $C_{ij} = \langle \varphi_i \varphi_j \rangle$ (with $C_{ii} = \langle \varphi_i^2 \rangle = \langle \varphi_i \rangle = \langle \varphi \rangle$).

Eq. (15) can be used as a starting point to gain deeper physical insight into the behavior of the standard deviation $\sigma$ than the previous result obtained directly from the grand partition function because it becomes possible to express $\sigma$ as an explicit function of $\langle \varphi \rangle$ and the correlation length $\xi$, see Eq. (35).

For $J > 0$, the standard deviation is bounded by the zero and infinite cooperativity values, which can be calculated from the limiting forms of $C_{ij}$. In the absence of cooperativity ($J = 0$), $C_{ij}$ factorizes for $i \neq j$: $C_{ij} \to \langle \varphi_i \rangle \langle \varphi_j \rangle = \langle \varphi \rangle^2$. In this $J = 0$ limit, $\langle \varphi^2 \rangle$ can then be obtained directly from Eq. (15) and $\sigma$ simplifies to $\sigma_0$ (Eq. (9)). In the other limit of infinite cooperativity ($J \to \infty$) and strong correlations, $\xi$ as a function of $\langle \varphi \rangle$ diverges at HF and away from HF for not too large $L$ remains larger than $L/2$ (see Fig. 15). From Eq. (29) we see that for $L$ finite in this limit $C_{ij} \to \langle \varphi_i^2 \rangle = \langle \varphi \rangle$, and $\sigma$ saturates at $\sigma_\infty$ (Eq. (12)).

### 2.3.4 Fluctuations in the strong anti-cooperativity limit

By comparing the $L = 13$ and $14$ cases (Figs. 3a and 3b), we observe that for small enough anti-cooperative systems, there is a clear difference between even and odd system sizes. This difference clearly manifests itself near HF, and, as we will see below, this difference arises because of frustration for odd $L$.

For negative values of $J$ (anti-cooperativity) and $L$ even the curves inflect and tend for $J \ll -1$ towards a characteristic limiting shape with global maxima near $\langle \varphi \rangle = 0.2$ and $0.8$ and a zero at half-filling (see Fig. 9). A restricted grand partition function approach (presented in Section 4.2), which retains only non-overlapping particle-hole pairs, leads to a simple but accurate approximation in the limit $J \to -\infty$ for $L$ even (see Fig. 9):

$$ \sigma_{-\infty}^{\text{even}} = L^{-1/2} \sqrt{\langle \varphi \rangle - 2\langle \varphi \rangle^2} \quad \text{for } 0 \le \langle \varphi \rangle \le 1/2. \qquad (16) $$

The corresponding standard deviation for $1/2 \le \langle \varphi \rangle \le 1$ can be obtained by exploiting particle-hole symmetry: $\sigma(\langle \varphi \rangle) = \sigma(1 - \langle \varphi \rangle)$.

In this limit of strong anti-cooperativity ($J \to -\infty$), the correlation length becomes complex, $1/\xi \approx i\pi + 2 \cosh(X)e^{-|J|/2}$. Hence, at HF, $1/\xi \approx i\pi$, leading directly to an oscillating two-point correlation function that describes a sequence of non-overlapping particle-hole pairs (anti-ferromagnetic order in Ising spin language):

$$ C_{i,i+r}^{\text{HF}} \approx \frac{1 + \cos(\pi r)}{4} \quad \text{for } J \to -\infty, \qquad (17) $$

where we have taken the real part of $C_{i,i+r}$ and assumed that $r/L \ll 1$. $C_{i,i+r}^{\text{HF}}$ alternates between 0 and $1/2$, which reflects alternating perfect anti-correlations (nearest neighbors) and correlations (next nearest neighbors) at HF in this limit of strong anti-cooperativity.

## 2.4 Probability distribution function (PDF) of the occupancy and its dependency on $\langle \varphi \rangle$ and $J$

We can extract the equilibrium probability distribution function (PDF) of the occupancy (corresponding to the effective Hamiltonian given in Eq. (1)) for fixed $J$ and $\mu$, $P(N; J, \mu)$ through exact enumeration and study how it depends on the interaction potential $J$. To study $P(N; J, \langle \varphi \rangle)$ for fixed $J$ and $\langle \varphi \rangle$, the experimentally accessible quantity, we varied $J$ but kept the mean occupancy fixed (by adjusting $\mu$ to any given choice of $J$).

Taking $L = 13$, Fig. 4 shows how the PDF behaves for three different positive values of the interaction potential and three characteristic values of the mean relative occupancy: $\langle \varphi \rangle = 1/3$, $1/2$, and $2/3$, respectively. The correlation length increases with $J$, and so does the standard deviation; consequently, the PDF profile broadens and flattens. On the one hand, finite-size effects set in when the tails of the PDF touch the system boundaries. On the other hand,
<figure><img src="image_6.png" /><figcaption>Figure 4: Discrete equilibrium probability distribution of the number of bound particles, P ( N ), for different values of the (dimensionless) interaction potential: J = 0 (orange), 3 (blue), and 5 (green). Histograms are obtained from an exact enumeration using the 1D Hamiltonian from Eq. (1) for a system of size L = 13 for the average occupancies 〈 N 〉 = 〈 ϕ 〉 L = 4 . 33, 6 . 5, and 8 . 67 ( 〈 ϕ 〉 = 1 / 3, 1 / 2, and 2 / 3), respectively.</figcaption></figure>

when the typical size of a highly correlated particle domain, ≈ 2 ξ , becomes comparable to the substrate size, L , finite-size effects dominate , and the PDF saturates to its strong interaction bimodal form.

2.4.1 PDF in the zero-cooperativity limit: the Hill-Langmuir case

In the absence of cooperativity, we expect N to be normally distributed with standard deviation Lσ$_{0}$ , 1 provided the system is large enough, and the position of the mean is far enough from the system boundaries: P$_{0}$ ( N ; 〈 ϕ 〉 ) ≡ P ( N ; 0 , 〈 ϕ 〉 ) ≈ N ( 〈 N 〉 , Lσ$_{0}$ ) for L ≫ 1. The half-filling case with J = 0 in Fig. 4b approximates this Gaussian behavior well because the average value, 〈 N 〉 , is far enough from the boundaries that the tails of the distribution do not reach the boundary values, 0 and L (the distribution decays fast enough to the left and right of the average value). More quantitatively, for any value of J , finite-size effects set in when 〈 N 〉 < 2 Lσ or L − 〈 N 〉 < 2 Lσ (i.e., when the tails of the PDF touch the system boundaries).

2.4.2 PDF in the strong cooperativity limit

The PDF broadens and flattens for intermediate positive values of the interaction potential. For higher values, e. g. J ≥ 5, the PDF saturates by accumulating at the boundaries at

$$
P_{\infty}(N; \langle \varphi \rangle) = (1 - \langle \varphi \rangle) \delta_{N,0} + \langle \varphi \rangle \delta_{N,L} \quad (18)
$$

(here, δ$_{ij}$ is the Kronecker delta) and the system becomes well described by an effective two-state (bimodal) system (this is the limit where there exists only one effective block domain). This result is intuitively plausible and can be obtained by retaining in the full partition sum only the empty and full states (see Section 4.2 for more details). In this limit, finite-size effects are always dominant, and the standard deviation for the occupancy N can easily be calculated directly from the PDF. As expected, it saturates at Lσ$_{∞}$ (= L/ 2 at half-filling). Owing to the simple form of this limiting PDF, it is easy to calculate all moments of N , leading to 〈 N $^{m}$〉 ∞ = L $^{m}$〈 ϕ 〉 ∞ (or 〈 ϕ $^{m}$〉 ∞ = 〈 ϕ 〉 ∞ ) and therefore explicit expressions for low order standardized moments, such as the skewness and kurtosis, see Section 4.2. This simplification occurs because, with the correlation length, the energy cost needed to create domain walls also increases with J . Hence, a strongly correlated system favors microscopic configurations that minimize the number of domain walls, selecting the mesostates corresponding to either an empty or a fully occupied

<footnote>$^{1}$Strictly speaking, the variable N being discrete, the probability follows a binomial distribution. But the latter is well approximated by the normal distribution in the large L limit.</footnote>
substrate. This tendency to select extremal occupancies is a characteristic signature of the presence of strong cooperativity effects in small-size systems. Furthermore, for a parameter choice corresponding to an occupancy expectation value different from half-filling, the distribution is skewed, as in Fig. 4a and 4c, with an asymmetry that increases with the value of the interaction potential.

### 2.4.3 PDF in the strong anti-cooperativity limit

The situation is more complicated for strong anti-cooperativity ($J \ll -1$) because the system cannot, in general, be reduced to a simple one or two-state system, except near zero and half-filling. Furthermore, near HF, the reduction depends on whether $L$ is even or odd (see Fig. 5). In the following, we will focus on the left half of the plots in Fig. 5, i. e. on occupancies from $N = 0$ to $L/2$, the right half being a reflection of the left through HF owing to the particle-hole symmetry of the system.

For $L$ even, the PDF has a single peak at $N = 0$ at zero filling and at $N = L/2$ at half-filling, leading to vanishing standard deviations in these two limiting cases (see Fig. 3b). At HF, there is only one allowed state, consisting of $L/2$ non-overlapping particle-hole pairs [$(\varphi_i, \varphi_{i+1}) = (1,0)$] (with a two-fold degeneracy, because the particles can all be on either the even or odd sites).

For $L$ odd, the PDF still has a single peak at $N = 0$ at zero filling, but since no single microstate corresponds to HF, there are in this case two peaks with equal (50%) weight at $N = (L-1)/2$ and $(L+1)/2$, leading to a non-vanishing standard deviation, $\sigma_{-\infty}^{\text{HF}} = 1/(2L)$, (see PDF for $L = 13$ in Fig. 5). For $L$ odd, particle-hole pairs cannot cover the whole system, and defects must appear to fulfill the HF constraint, either an extra hole ($N = (L-1)/2$) or an extra particle ($N = (L+1)/2$).

### 2.4.4 PDF near half-filling: summary

We can summarize the situation as follows. At (or close to) HF, the four different studied cases give rise to very different results for the fluctuations described by the standard deviation and can be considered as clearly distinguishable signatures of different types of particle-particle correlations in small-size systems. At half-filling, the standard deviation takes on the following values in decreasing order: (i) for strong cooperativity ($J \to \infty$), $\sigma = 1/2$; (ii) for no cooperativity ($J = 0$), $\sigma = 1/(2\sqrt{L})$; (iii) for strong anti-cooperativity ($J \to -\infty$) and $L$ odd, $\sigma = 1/(2L)$; and (iv) for strong anti-cooperativity ($J \to -\infty$) and $L$ even, $\sigma = 0$.

In contradistinction to the strong cooperativity limit where the two-state approximation is valid over the whole range of relative occupancy (from zero to full-filling) for both $L$ even and odd, the strong anti-cooperativity limit is more complicated because the one or two-state approximation is only valid for zero and half-filling. Although more than two occupancy states are involved between these limits, it is clear that strong anti-cooperativity favors the formation of particle-hole pairs. Therefore for even $L$, the system can be approximated by a system of non-overlapping particle-hole pairs. For $L$ odd, two extra defect states must be included (an extra particle or extra hole). This restricted grand partition function approach, where only the relevant states participating in the strong cooperativity and strong anti-cooperativity limits are retained, is developed in Section 4.2 of the Methods.

## 2.5 Application to the bacterial flagellar motor (BFM)

In what follows, we aim to determine whether cooperativity between BFM stator units, mediated by a (dimensionless) interaction potential $J$, plays a role in their dynamical assembly at the periphery of the rotor.

Like most microorganisms, bacteria live in fluid environments with low Reynolds numbers, making them experience a viscous force much larger than the inertial ones [20]. They have
<figure><img src="image_7.png" /><figcaption>Figure 5: Discrete equilibrium probability distribution of the number of bound particles, $P(N)$, for $J = -10$ and three different mean occupancies: $\langle\phi\rangle = 1/3, 1/2$, and $2/3$, respectively. Histograms show the result of exact enumeration based on the 1D Hamiltonian from Eq. (1) for a system of size $L = 13$ (a) and $L = 14$ (b).</figcaption></figure>

subsequently evolved a variety of compatible motility mechanisms that have been widely studied [21–26]. One type of motility involves rotating one or several flagella that propel bacteria through aqueous media [27–35]. The BFM is a transmembrane macromolecular complex that consumes the electrochemical potential across the inner bacterial cell membrane to generate torque and set the flagellum in rotary movement [8, 9, 12, 31, 36, 37]. One of the most compelling properties of the BFM is its ability to change both its conformation and stoichiometry depending on the external medium, allowing it to change the direction of rotation and the magnitude of torque produced [11, 13, 38, 39]. Its mechanisms of adaptation to external stimuli have been widely studied, becoming a model molecular machine to investigate properties such as mechanosensitivity [11, 14, 39–42], chemotaxis [43–45] and dynamic subunit exchange [46].

Torque is generated by inner membrane ion channel complexes called stator units which dynamically bind and unbind to the peptidoglycan (cell wall) at the periphery of the rotor. When unbound, they are inactive and passively diffuse in the inner membrane. In their bound state, anchored to the peptidoglycan, the ion channels are activated, and, through a mechanism not yet fully understood, apply torque to the rotor [47–54]. Precise measurements of the temporal evolution of the angular velocity of the motor of *E. coli*, a direct proxy of the number of bound stator units, have shown that the system can recruit up to $L \approx 13$ stator units and that the system is mechanosensitive in that stator unit stoichiometry scales with the external torque induced by the viscous drag acting on the flagellum turning inside an aqueous medium [34].

Much can be learned about the dynamics of the BFM using bead assay measurements. In our experiments performed on *E. coli*, we attach a microparticle (bead) to the ‘hook’ (the extracellular portion that joins the motor to the flagellum) via a flagellum stub. By tracking the bead’s off-axis rotation, we can calculate the angular velocity, $\omega$, and the torque produced, $\tau$ (from the relation $\tau = \gamma\omega$, where $\gamma$ is the drag coefficient that increases with the bead’s diameter), both of which are a direct proxy for the number of bound stator units. We can (indirectly) control the mean number of bound stator units at steady state, $\langle N \rangle$, by varying the beads’ size and hence the viscous load, because the binding of stators to the BFM is mechanosensitive (see, e.g., [40], and references therein). We can thus measure the temporal evolution of the number of bound stator units on individual motors, as well as the fluctuations around mean occupancy (for more information on the experimental setup, see references [40, 55]).

We describe the mechanosensitive binding and unbinding of stator units in the stationary angular velocity regime of the BFM in terms of our adsorption model (see Eq. (1)), the rotor being a small-size substrate with periodic boundary conditions onto which up to $L = 13$ stator units can bind at fixed equally spaced positions. Our working hypothesis, in line with previous studies, is that as long as we only focus on the occupancy of the stator units in the (non-
<figure><img src="image_8.png" /><figcaption>Figure 6: Parametric plot of the standard deviation of the relative occupancy, $\sigma$, versus the mean occupancy, $\langle\varphi\rangle$, comparing the experimental values with the theory; here, $L = 13$. The experimental data corresponds to five different applied viscous loads with three different micro-particles (beads) whose diameter is indicated in the legend. The experimental data fall between the values of $J = 0.5$ and $J = 2$ of the theoretical curves (aside from one outlier). From this data, we can estimate a characteristic experimental error of $\Delta\sigma \approx 0.02$.</figcaption></figure>

equilibrium) motor stationary state this quantity fluctuates around a fixed mean value as if the stator subsystem were effectively at equilibrium. In this picture, we account for the presence of the unbound (inactive) stator units diffusing freely in the inner membrane by imposing an external reservoir chemical potential, $\mu_r$, which is taken to be constant, as we expect depletion effects to be negligible. We incorporate the mechanosensitivity into the model in an average way by assuming that the stator unit binding energy, $\epsilon$, depends on the viscous load (in our case dependent on the size of the bead and the viscosity of the surrounding medium). Load-induced changes in $\epsilon$ will naturally lead to load-dependent average occupancies and fluctuations, as observed experimentally. Furthermore, we assume that the interaction parameter $J$ remains fixed (independent of the load) and check a *posteriori* if this assumption is consistent with the data.

As explained previously, we attempt to use the fluctuations in the average number of bound stator units in the stationary angular velocity regime to determine whether or not cooperativity is at play in the BFM. Fig. 6 shows a comparison between the theoretical standard deviation curves, already presented in Fig. 3a, and the experimental standard deviations from five different applied viscous loads, corresponding to three different beads with different diameters that are indicated in the legend. One can see that the experimental data fall between the values of $J = 0.5$ and 2, leading to the conclusion that (i) a constant (load-independent) cooperativity parameter ($J$) is a reasonable working hypothesis, (ii) a moderate level of cooperativity in the system, $J \approx 1.21 \pm 0.22$, obtained by fitting the model to the experimental data, is consistent with the experimental observations, see Fig. 17 (0.22 is the standard error of the nonlinear fit); and (iii) the estimated value of $J$ is coherent with what is expected for typical biological systems exhibiting moderate cooperativity.

To delve into further detail, we can also compare the probability distributions of the occupancy shown in Section 2.4 with our preliminary experimental ones. Fig. 7 shows a comparison between the probability distribution from experiments (black dots) and the probability distributions obtained by exact enumeration for different values of $J$. Even though, in this analysis of preliminary data, we don't observe a precise fit with any of the chosen values, it is clear that the PDFs for the range of $J$ estimated from the standard deviation are coherent with our experimental results. This result suggests that the interaction potential has a non-zero positive value qualitatively in accordance with the value range deduced from the above fluctuation analysis. Moreover, as further explained in the Methods section, we observe theoretically that the system becomes bimodal (with only the completely full and empty states having significant
<figure><img src="image_9.png" /><figcaption>Figure 7: Discrete equilibrium probability distribution of bound stator units, $P(N)$, for the experiments (red with black dots) and different values of the interaction potential, $J = 0$ (blue), 1 (orange), 2 (green), obtained from exact enumeration of states from the effective Hamiltonian; here, $L = 13$. The mean equilibrium occupancy of each distribution is indicated on top of the three plots, which corresponds to a given external load imposed by beads of different diameters indicated in the legend.</figcaption></figure>

probabilities) for slightly higher values of the interaction potential, $J \gtrsim 5$, behavior that is not seen in the experimental PDF data. We see from Section 2.4 (Fig. 4) that for the chosen system size ($L = 13$), $J \approx 3$ is a threshold value marking a transition from low cooperativity PDFs having a maximum centered on the average occupancy to high cooperativity PDFs exhibiting a bimodal form with local maxima at zero and full filling. In the latter case, microstates with fillings close to the average occupancy have low probabilities, and the motor would undergo discontinuous jerky motion. Such a bacteria would fluctuate between an immobile state (zero stator unit occupancy) and maximum speed (full stator unit occupancy). Knowing that the BFM is a highly and smoothly adaptive molecular machine, such high values of cooperativity would, therefore, not be expected since such a motor would not be able to adapt smoothly to environmental variations.

## 3 Discussion and Conclusion

We have focused on cooperative processes involving the adsorption of ligands onto a substrate disposing of a limited number of binding sites. By presenting a general method (based on a 1D periodic lattice gas) for recognizing and assessing characteristic signatures of cooperativity or anti-cooperativity in the stochastic occupancy fluctuations, we propose both a criterion to determine whether any given adsorption system exhibits cooperative or anti-cooperative behavior and a method to quantify the amplitude of the ligand-ligand interaction potential.

In the process, knowing that in the thermodynamic limit relative occupancy fluctuations (or standard deviation) vanish, we have addressed the following essential questions: (i) what is a sufficiently “small” system for studying fluctuations? (ii) what model parameter values allow a system to smoothly “adapt” to external conditions?

We compared the theoretical results for both the standard deviation and the probability distribution function of stator unit occupancy with experimental data obtained for the BFM. We concluded that a moderate value of cooperativity, ≈ 1 $k_B T$ (i.e., $J \approx 1$ with a standard error of about 20%), for the short-range (SR) model is not only coherent with the experimental data but also expected from the characteristic smooth adaptability of the motor to changing external loads [for the infinite-range (IR) model, using the mapping (72) derived in Section 4.3, the corresponding value is $J_{IR} = 2J_{SR}/(L-1) \approx 0.20$]. For the characteristic size of the BFM (≈ 10 binding sites), slightly higher values of cooperativity would lead to a motor that is bimodal in the occupancy and, therefore, would exhibit switch-like behavior for the produced torque not compatible with the required motor characteristics. The highly cooperative bimodal behavior
found at strong coupling does, however, appear to describe the rapid stochastic switching in BFM rotational direction [11, 17, 56].

Within the framework that we have developed, we are now in a position to propose, using the BFM as an illustrative example, a general principle of motor *adaptability* depending on whether the motor under investigation should respond smoothly to external stimuli or behave like a two-state switch. As stated earlier, we suppose that for the BFM the binding energy $\epsilon$ depends on the load (and therefore the bead size). A change in load would therefore lead directly to a change in the effective chemical potential $\mu$ and therefore a modification of the average occupancy and standard deviation (fluctuations).

For the short-range model, these modifications to $\mu$ (extracted from the experimental data) are shown in Fig. 8 on the occupancy and standard deviation contour plots for each of the three studied values of $J$ (for $L = 13$). By positioning the effective chemical potential $\mu$ window in this way we observe that for $1 < J < 2$ it is situated in a *sweet spot* suitable for a motor that responds smoothly to environmental changes with the ability to cover a wide range of occupancies spanning half-filling while minimizing the amplitude of the fluctuations.

For the BFM the occupancy is directly related to the motor speed and therefore, near this *sweet spot*, the BFM can smoothly adjust its speed in response to external stimuli with a minimum of fluctuations. In retrospect the bead sizes used in the experiments were clearly chosen to *see* an effect because a much higher load would have forced the system into nearly full filling (exactly the case in the stall experiments reported in [55]); a much lighter load would have pushed the system to nearly zero filling (similar to what was done using another technique in the resurrection experiments). The highly cooperative switch-like behavior observed for the BFM rotational direction [11, 56] would place this system in the strong coupling regime on the right-hand side of the contour plots ($J > 4$).

For the motor to respond smoothly to external stimuli *and* cover a wide range of occupancy with minimum of fluctuations we therefore see that moderate positive values of cooperativity $1 < J < 2$ are optimal given the $\mu$ window imposed by the system characteristics. Although the occupancy range would increase at higher values of $J$, this positive effect would be counterbalanced by a strong increase in fluctuations because the system would be pushed into the switch-like operation regime. On the other hand, for lower and even negative values of $J$, the fluctuations would (favorably) be diminished in amplitude, but at the cost of severely restricting the accessible range of occupancy (and therefore motor speed for the BFM).

By examining Fig. 8, we observe that if a two-state switch-like motor operation were sought after in order to cover a wider range of average relative occupancy, then higher values of $J$ would be the best choice, the price to pay would be a strong increase in fluctuations, provided that the observation time window be wide enough. If a motor with weak fluctuations and a restricted range of occupancies were sought after then negative values of $J$ (anti-cooperativity) would be the best choice. In fact, this type of motor would be relatively immune to changes in external stimuli and therefore exhibit a relatively constant speed (provided the speed were still directly related to occupancy, as for the BFM).

We have not addressed here the biochemical origins of stator-stator interactions, although one can imagine that stator units interact at short range much like proteins, either directly or through allosteric pathways, and that long-range interactions could also be due to allosteric effects. We plan to address this and other open questions in future work. In the present context the major important open question in modeling the BFM concerns how to integrate cooperativity into kinetic models that already allow one to account for the relaxation time asymmetry between stall and resurrection [55]. The goal is to arrive at a model that can describe in a unified way both BFM stator number fluctuations and relaxation time asymmetries.
<figure><img src="image_10.png" /><figcaption>Figure 8: Values of the effective chemical potential $\mu$ extracted from the experimental data for different values of the interaction potential $J$ ($L = 13$). We situate them on the contour plots of $\langle\varphi\rangle$ and $\sigma$ in order to find the interval of chemical potential for which the BFM would work at fixed $J$.</figcaption></figure>

# 4 Methods

## 4.1 1D short-range lattice gas (SRLG)

Analytical solutions for the equilibrium state can be obtained using the transfer matrix formalism [18, 57]. According to it, the grand partition function of the system can be written as:

$$ \Xi = \sum_{\{\varphi_i\}} \exp[-\beta H(\varphi)] = \sum_{\{\varphi_i\}} \prod_{i=1}^{L} \langle\varphi_i|T|\varphi_{i+1}\rangle, \quad (19) $$

where $T$ is the transfer matrix, which for the system we are working on, is given by

$$ T = \begin{pmatrix} 1 & e^{\frac{\mu}{2}} \\ e^{\frac{\mu}{2}} & e^{(J+\mu)} \end{pmatrix}. \quad (20) $$

In the case of periodic boundary conditions, equation (19) can be simplified to

$$ \Xi = \mathrm{Tr} \, T^L = \lambda_+^L + \lambda_-^L \quad (21) $$

where

$$ \lambda_{\pm} = e^X \left( \cosh X \pm \sqrt{\sinh^2 X + e^{-J}} \right) \quad (22) $$

are the eigenvalues of the transfer matrix and $X = \frac{1}{2}(J + \mu)$. From this result (21) for the grand partition function, one can obtain the mean relative occupancy, $\langle\varphi\rangle$, and the mean square relative occupancy, $\langle\varphi^2\rangle$, as functions of $\mu$ and $J$ by employing known relations between thermodynamic quantities (i.e., by taking derivatives of $\Xi$ with respect to the effective chemical potential $\mu$). The standard deviation $\sigma$ can then be derived as a function of $\mu$ and $J$ from these two averages.

The transfer matrix formalism can also be used to calculate the occupancy-occupancy correlation function, which provides another route to $\langle\varphi^2\rangle$ (and therefore $\sigma$), but now as an explicit function of $\langle\varphi\rangle$ and the correlation length, $\xi \equiv 1/\ln(\lambda_+/\lambda_-)$ (already presented in Eq. 3). This second form for $\sigma$ lends itself to a clearer physical picture of how the standard deviation interpolates between the zero and infinite cooperativity limits.

According to the grand canonical formalism, the mean relative occupancy of a system is given by the expression

$$ \langle\varphi\rangle = \frac{1}{L} \frac{\partial \Xi}{\partial \mu} = - \frac{1}{L^2} \frac{\partial \Omega}{\partial \mu}, \quad (23) $$
where $\Omega = -L \ln \Xi$ is the grand potential corresponding to the grand partition function. Using the chain rule, we can find that

$$ \langle \varphi \rangle = \frac{1}{\Xi} \left( \lambda_{+}^{(L-1)} \frac{\partial \lambda_{+}}{\partial X} \frac{\partial X}{\partial \mu} + \lambda_{-}^{(L-1)} \frac{\partial \lambda_{-}}{\partial X} \frac{\partial X}{\partial \mu} \right). \quad (24) $$

Replacing

$$ \frac{\partial X}{\partial \mu} = \frac{1}{2} \quad \text{and} \quad \frac{\partial \lambda_{\pm}}{\partial X} = \lambda_{\pm} \left( 1 \pm \frac{\sinh X}{\sqrt{\sinh^2 X + e^{-J}}} \right) \quad (25) $$

in the previous equation, we obtain Eq. 2.

Similarly, we can determine the mean square relative occupancy by taking another derivative:

$$ \langle \varphi^2 \rangle = \frac{1}{L^2} \frac{1}{\Xi} \frac{\partial^2 \Xi}{\partial \mu^2}. \quad (26) $$

By carrying out this calculation, we obtain:

$$ \langle \varphi^2 \rangle = \frac{1}{4} \left[ 1 + \frac{\sinh^2 X}{\sinh^2 X + e^{-J}} + \tanh \left( \frac{L}{2\xi} \right) \left( \frac{2 \sinh X}{\sqrt{\sinh^2 X + e^{-J}}} + \frac{1}{L} \frac{e^{-J} \cosh X}{(\sinh^2 X + e^{-J})^{3/2}} \right) \right]. \quad (27) $$

Denoting $\xi$ as the correlation length makes physical sense because correlations decay exponentially on a scale $\xi$, as can be seen by inspecting the correlation function obtained using the transfer matrix formalism:

$$ \langle \varphi_i \varphi_{i+r} \rangle = \frac{1}{K} \left[ \lambda_{+}^{L} (1 + \lambda_{+}^{-r} \lambda_{-}^{r} v_{-}^{2}) + \lambda_{-}^{L} v_{-}^{2} (v_{-}^{2} + \lambda_{+}^{r} \lambda_{-}^{-r}) \right], \quad (28) $$

where $v_{-} = e^{\mu/2}(\lambda_{-} - 1)$ and $K = (\lambda_{+}^{L} + \lambda_{-}^{L})(1 + v_{-}^{2})^{2}$. By using the definition of $\xi$ presented above, the products of $\lambda_{+}^{r} \lambda_{-}^{-r}$ and $\lambda_{+}^{-r} \lambda_{-}^{r}$ can be rewritten as $e^{\pm r/\xi}$. This result for the correlation function, Eq. 28, can be put into a more physically transparent form,

$$ \langle \varphi_i \varphi_{i+r} \rangle = \langle \varphi \rangle + \frac{\langle \varphi \rangle_{\text{TL}} - \langle \varphi \rangle_{\text{TL}}^2}{1 + e^{-L/\xi}} \left[ e^{-r/\xi} + e^{-(L-r)/\xi} - 1 - e^{-L/\xi} \right]. \quad (29) $$

This result, which illustrates how the correlation function can be written explicitly in terms of $\langle \varphi \rangle$, $\xi$, and $L$, allows the exponential decay on a scale $\xi$ to be made clear. Here $\langle \varphi \rangle_{\text{TL}}$ is the average relative occupancy in the thermodynamic limit (TL), defined to be an explicit function of $\langle \varphi \rangle$ and $\xi$:

$$ \langle \varphi \rangle_{\text{TL}} \equiv \frac{1}{2} \left[ 1 + \coth \left( \frac{L}{2\xi} \right) (2\langle \varphi \rangle - 1) \right]. \quad (30) $$

This result for the correlation function, Eq. (29), is clearly periodic and tends to the correct TL when $L \to \infty$:

$$ \langle \varphi_i \varphi_{i+r} \rangle_{\text{TL}} = \langle \varphi \rangle_{\text{TL}}^2 + \left( \langle \varphi \rangle_{\text{TL}} - \langle \varphi \rangle_{\text{TL}}^2 \right) e^{-r/\xi}. \quad (31) $$

By performing the sum depicted in Eq. 15 using Eq. 28 we obtain for the mean square relative occupancy an expression that explicitly depends on the eigenvalues and the correlation length as follows:

$$ \langle \varphi^2 \rangle = \frac{1}{K} \left[ \lambda_{+}^{L} \left( 1 + \frac{v_{-}^{2}}{L} \frac{1 - e^{-L/\xi}}{1 - e^{-1/\xi}} \right) + \lambda_{-}^{L} v_{-}^{2} \left( v_{-}^{2} + \frac{1}{L} \frac{1 - e^{L/\xi}}{1 - e^{1/\xi}} \right) \right]. \quad (32) $$

It can be shown that Eq. 32 is equivalent to Eq. 27, which provides a check on the thermodynamic result. In the calculation leading to Eq. 27, the correlation length appears from the fraction $(\lambda_{+}^{L} - \lambda_{-}^{L})/(\lambda_{+}^{L} + \lambda_{-}^{L})$ which can be rewritten as $\tanh(L/(2\xi))$.
With both Eq. 2 and Eq. 27, we can obtain an analytical expression for the standard deviation from the expression $\sigma = \sqrt{\langle\varphi^2\rangle - \langle\varphi\rangle^2}$. The result obtained in this way is equivalent to the one obtained more efficiently in the main text using the thermodynamic relation $\sigma = L^{-1/2}\sqrt{\partial_\mu\langle\varphi\rangle}$.

To make explicit the $\xi$ dependence we now calculate $\sigma$ using the connected correlation function, given by $\langle\varphi_i\varphi_j\rangle - \langle\varphi\rangle^2$, in the double sum leading to $\sigma^2$:

$$ \sigma^2 = \langle\varphi^2\rangle - \langle\varphi\rangle^2 = \frac{1}{L^2} \sum_{i=1}^{L} \sum_{j=1}^{L} [\langle\varphi_i\varphi_j\rangle - \langle\varphi\rangle^2] \quad (33) $$

$$ = \frac{\langle\varphi\rangle}{L} - \langle\varphi\rangle^2 + \frac{2}{L} \sum_{r=1}^{L-1} \left(1 - \frac{r}{L}\right) \langle\varphi_1\varphi_{1+r}\rangle. \quad (34) $$

The sum in 34 can be evaluated using Eq. 29. to obtain a physically transparent expression for the variance, $\sigma^2$, explicitly in terms of $\langle\varphi\rangle$, $\xi$, and $L$:

$$ \sigma^2 = \sigma_{\infty}^2 + (\langle\varphi\rangle_{\text{TL}} - \langle\varphi\rangle_{\text{TL}}^2) F(\xi; L), \quad (35) $$

where $\langle\varphi\rangle_{\text{TL}}$ is given by Eq. 30 and

$$ F(\xi; L) \equiv \frac{(L+1)\left(e^{-\frac{1}{\xi}} - e^{-\frac{L}{\xi}}\right) + (L-1)\left(e^{-\frac{(L+1)}{\xi}} - 1\right)}{L\left(1 - e^{-\frac{1}{\xi}}\right)\left(1 + e^{-\frac{L}{\xi}}\right)}. \quad (36) $$

This form for the variance allows the weak and strong cooperativity limits to be clearly identified, since $F \to 0$, when $J \to \infty$ and $F \to -1 + 1/L$, when $J \to 0$ (in the latter limit, $\langle\varphi\rangle_{\text{TL}} \to \langle\varphi\rangle$). More generally, for fixed $J$, $\xi$ becomes an implicit function of $\langle\varphi\rangle$ and can be plotted parametrically, as in Fig. 15 for infinite $J$.

The thermodynamic limit corresponds to $L \to \infty$, for which $\tanh(L/2\xi) \to 1$, therefore the average relative occupancy simplifies to

$$ \langle\varphi\rangle_{\text{TL}} = \frac{1}{2} \left( 1 + \frac{\sinh X}{\sqrt{\sinh^2 X + e^{-J}}} \right), \quad (37) $$

and $\langle\varphi^2\rangle_{\text{TL}} = \langle\varphi\rangle_{\text{TL}}^2$. Hence the standard deviation in the thermodynamic limit, $\sigma_{\text{TL}}$, is zero and self-averaging holds.

## 4.2 Probability distribution function of the occupancy

For high positive values of $J$, e. g. $J \ge 5$, the PDF saturates by accumulating at the boundaries to

$$ P_{\infty}(N; \langle\varphi\rangle) = (1 - \langle\varphi\rangle)\delta_{N,0} + \langle\varphi\rangle\delta_{N,L} \quad (38) $$

and the system becomes well described by an effective two-state system. The standard deviation for the occupancy $N$ therefore saturates at $L\sigma_\infty (= L/2$ at half-filling), where we recall that $\sigma_\infty = \sqrt{\langle\varphi\rangle - \langle\varphi\rangle^2}$, and it becomes easy to calculate all moments of $N$: $\langle N^m \rangle = L^m \langle\varphi\rangle$ (or $\langle\varphi^m\rangle = \langle\varphi\rangle$). The moments are defined by

$$ \mu_n = \langle (\varphi - \langle\varphi\rangle)^n \rangle \quad (39) $$

and the standard moments by $\mu_n/\sigma^n$, where $\sigma = \sqrt{\mu_2}$ is the standard deviation. In the strong cooperativity limit, one can therefore find simple explicit expressions for low order standard moments, such as the skewness $\gamma = \mu_3/\sigma^3$,

$$ \gamma_{\infty} = \frac{\mu_{3,\infty}}{\sigma_{\infty}^3} = \frac{\langle\varphi\rangle - 3\langle\varphi\rangle^2 + 2\langle\varphi\rangle^3}{\sigma_{\infty}^3} \quad (40) $$
and kurtosis, $\kappa = \mu_4/\sigma^4$,

$$ \kappa_{\infty} = \frac{\mu_{4,\infty}}{\sigma_{\infty}^4} = \frac{\langle\varphi\rangle - 4\langle\varphi\rangle^2 + 6\langle\varphi\rangle^3 - 3\langle\varphi\rangle^4}{\sigma_{\infty}^4}. \qquad (41) $$

These results for the effective two-state system are very different from those predicted for a Gaussian PDF.

For strong anti-cooperativity ($J \ll -1$), the system near zero and half-filling can be reduced to a simple one or two-state system. Near HF, the reduction depends on whether $L$ is even or odd (see Fig. 5). The PDF always has a single peak at $N = 0$ at zero filling. For even $L$ in the strong anti-cooperativity limit, the PDF also has a single peak at $N = L/2$ at half-filling,

$$ P_{-\infty}(N; 1/2) = \delta_{N,L/2} \quad (L \text{ even}), \qquad (42) $$

(see Fig. 5) leading to a vanishing standard deviation. For $L$ odd at half-filling, on the other hand, there are two peaks with equal (50%) weight at $N = (L - 1)/2$ and $(L + 1)/2$,

$$ P_{-\infty}(N; 1/2) = \frac{1}{2}\delta_{N,(L-1)/2} + \frac{1}{2}\delta_{N,(L+1)/2} \quad (L \text{ odd}), \qquad (43) $$

leading to a non-vanishing standard deviation of $1/2$ for the occupancy $N$. Unlike for $L$ even, for $L$ odd non-overlapping particle-hole pairs cannot cover the whole system, and a defect (non-particle-hole pair) must appear, either an extra hole [$N = (L - 1)/2$] or an extra particle [$N = (L + 1)/2$], to arrive at HF.

The restricted grand partition function approach consists in keeping in the sum over all possible microstates $\{\varphi_i\}$ only those states that survive in the studied limit. Before presenting this approach, we first recall the exact calculation in the absence of cooperativity ($J = 0$).

For $J = 0$ the partition sum can be organized into a sum over states with a fixed occupancy $N$ with the Boltzmann factor $\exp[-\beta H_0(N)] = e^{\mu N}$ multiplied by a multiplicity (or binomial coefficient) $C_N^L = L! / [N! (L - N)!]$ (related to the configurational entropy) that gives the number of microstates $\varphi$ consistent with occupancy $N$:

$$ \Xi_0 = \sum_{\{\varphi_i\}} e^{-\beta H_0 \{\varphi_i\}} = \sum_{N=0}^{L} C_N^L e^{\mu N} = (1 + e^{\mu})^L. \qquad (44) $$

From $\Xi_0$ we immediately obtain the expected $J = 0$ results:

$$ \langle\varphi\rangle_0 = (L\Xi_0)^{-1}\partial_{\mu}\Xi_0 = 1/(1+e^{-\mu}) \qquad (45) $$

and

$$ \sigma_0 = \sqrt{L^{-1}\partial_{\mu}\langle\varphi\rangle_0} = L^{-1/2}\sqrt{\langle\varphi\rangle - \langle\varphi\rangle^2}. \qquad (46) $$

By introducing the usual (configurational) entropy of mixing,

$$ S_{\text{mix}} = -k_B L[\varphi \ln(\varphi) + (1-\varphi) \ln(1-\varphi)], \qquad (47) $$

the multiplicity $C_N^L$ can be accurately approximated (using Stirling's formula) for $L \gg 1$ and $0 < \varphi = N/L < 1$ by

$$ C_N^L \approx \frac{e^{Ls_{\text{mix}}(\varphi)}}{\sqrt{2\pi L \varphi (1-\varphi)}}, \qquad (48) $$

where

$$ s_{\text{mix}} \equiv S_{\text{mix}}/(Lk_B) = -[\varphi \ln(\varphi) + (1-\varphi) \ln(1-\varphi)] \qquad (49) $$

is the dimensionless entropy of mixing. The entropy of mixing is a non-monotonic (concave) function of $\varphi$ that goes to 0 for zero and full filling and reaches a maximum at HF.
In Eq. (48) for $L \gg 1$, the main variation comes from the exponential, and it is possible under certain conditions to replace $\varphi$ in the prefactor by $\langle \varphi \rangle$ without loss of accuracy. The PDF for $J = 0$ can then be accurately approximated by

$$P_0(N; \langle \varphi \rangle) \equiv \frac{C_N^L e^{\mu N}}{\Xi_0} \approx \frac{e^{-\beta \mathcal{F}_0(\varphi; \langle \varphi \rangle)}}{(2\pi)^{1/2} \sigma_0 L \Xi_0(\langle \varphi \rangle)}, \quad (50)$$

where

$$\beta \mathcal{F}_0(\varphi; \langle \varphi \rangle) = -L \left[ \varphi \ln \left( \frac{\langle \varphi \rangle}{1 - \langle \varphi \rangle} \right) + s_{\text{mix}}(\varphi) \right] \quad (51)$$

is an effective free energy for $\varphi$ and we have written

$$\Xi_0(\langle \varphi \rangle) = \frac{1}{(1 - \langle \varphi \rangle)^L}, \quad (52)$$

and

$$\mu(\langle \varphi \rangle) = \ln \left( \frac{\langle \varphi \rangle}{1 - \langle \varphi \rangle} \right) \quad (53)$$

as functions of $\langle \varphi \rangle$.

For $L \gg 1$, we can use the saddle point approximation to obtain a simple Gaussian form for $P_0(N; \langle \varphi \rangle)$: by expanding $\beta \mathcal{F}_0(\varphi; \langle \varphi \rangle)$ around $\varphi = \langle \varphi \rangle$ (the position of its minimum) to second order in $(\varphi - \langle \varphi \rangle)$ we find

$$P_0(N; \langle \varphi \rangle) \approx \frac{\exp \left[ - \frac{(\varphi - \langle \varphi \rangle)^2}{2\sigma_0^2} \right]}{(2\pi)^{1/2} \sigma_0 L}, \quad (54)$$

where the variance

$$\sigma_0^2 = -[L s''_{\text{mix}}(\langle \varphi \rangle)]^{-1} = (\langle \varphi \rangle - \langle \varphi \rangle^2)/L \quad (55)$$

is related to the second derivative of the entropy of mixing, or curvature, evaluated at $\langle \varphi \rangle$. Since $-s''_{\text{mix}}(\langle \varphi \rangle)$ has a minimum at half-filling, $\sigma_0$ has a maximum there.

The Gaussian approximation for $P_0(N; \langle \varphi \rangle)$ is (approximately) correctly normalized because the sum of the exact PDF over $N$ is equal to 1, which is equivalent to the Gaussian approximation times $L$ integrated over $\varphi$ from $-\infty$ to $+\infty$ about its maximum value at $\langle \varphi \rangle$. The Gaussian approximation is itself clearly only valid when the tails of the distribution are far enough from the extremal values ($\varphi = 0$ and 1), which is the case for the three $J = 0$ PDFs displayed in Fig. 4 (for $\langle \varphi \rangle = 0.31$, 0.5, and 0.62). For these three cases the Gaussian approximation is extremely accurate. More generally a PDF calculated using the full approximation for $C_N^L$, Eq. (48), should be accurate over the whole range of $\varphi$ except close to the extremal values.

In the strong cooperativity limit ($J \to +\infty$), only the empty and full states need be retained in the restricted grand partition function, $\Xi_{+\infty}$, because any state with $0 < N < L$ will necessarily have domain walls and these states will be suppressed via the Boltzmann factor by $e^{-J}$ (with respect to the fully occupied state) for each lost nearest-neighbor interaction. This leads to

$$\Xi_{+\infty}(X; \mu) = 1 + e^{2XL}, \quad (56)$$

where $X = (J + \mu)/2$ is considered to be kept fixed as $J \to +\infty$ and $\mu \to -\infty$. From $\Xi_{+\infty}(X; \mu)$ we obtain

$$\langle N \rangle_{+\infty} = \Xi_{+\infty}^{-1} \partial_{\mu} \Xi_{+\infty} = L e^{2XL} / \Xi_{+\infty}, \quad (57)$$

which allows us to read off the expected results,

$$P_{+\infty}(0; \langle \varphi \rangle) = \Xi_{+\infty}^{-1} = 1 - \langle \varphi \rangle \quad (58)$$

and

$$P_{+\infty}(L; \langle \varphi \rangle) = e^{2XL} \Xi_{+\infty}^{-1} = \langle \varphi \rangle, \quad (59)$$
directly from $\Xi_{+\infty}$ [within this approximation the PDF $P_{+\infty}(N; \langle\varphi\rangle)$ vanishes for all other values of $N$].

For strong anti-cooperativity ($J \ll -1$), we must treat the $L$ odd and even cases differently. In both these cases, the statistical physics is very different from the strong cooperativity ($J \gg 1$) limit where only two states need to be retained.

We start with the simpler $L$ even case for which the system can be approximated by non-overlapping, and therefore non-interacting particle-hole pairs (we choose the convention where the position of the particle-hole pair is determined by the position of the particle making up the particle-hole pair). The restricted grand partition function for $0 \le \langle\varphi\rangle \le 1/2$ in this particle-hole pair approximation is

$$ \Xi_{-\infty}^{\text{even}} = \sum_{N_d=0}^{L/2} D_{N_d}^{L} e^{\mu N_d}, \qquad (60) $$

where the upper limit on the sum of $L/2$ reflects the maximum number of non-overlapping pairs on a lattice of size $L$ (even) and the multiplicity $D_{N_d}^L$ counts the number of distinct ways of putting $N_d$ pairs on a periodic lattice of size $L$. Since calculating $D_{N_d}^L$ is a non-trivial combinatorial problem, we simplify matters by considering two disjoint lattices of size $L/2$ consisting of the odd sites for the first and the even sites for the second. If we neglect mixed particle-hole pair microstates, i.e., those that have pairs on both lattices, we obtain the following approximate restricted particle-hole pair grand partition function

$$ \Xi_d^{\text{even}} = \sum_{N=0}^{L/2} 2C_{N_d}^{L/2} e^{\mu N_d} = 2(1 + e^{\mu})^{L/2}, \qquad (61) $$

where $C_{N_d}^{L/2}$ is the usual binomial coefficient that counts the number of distinct ways of placing $N_d$ particles on a lattice of size $L/2$. In passing from Eq. (60) to Eq. (61), we have replaced $D_{N_d}^L$ by $2C_{N_d}^{L/2}$ and therefore simplified the problem to that of two independent non-interacting quasi-particle (i.e., particle-hole pair) systems [cf. Eq. (44)]. Following the above discussion for the $J=0$ case, we can immediately conclude that $\langle\varphi\rangle = \langle N_d/L \rangle \approx \langle\varphi\rangle_d^{\text{even}}$, where

$$ \langle\varphi\rangle_d^{\text{even}} = (L\Xi_d)^{-1}\partial_{\mu}\Xi_d^{\text{even}} = 1/[2(1+e^{-\mu})] \qquad (62) $$

and therefore $\sigma_{-\infty}^{\text{even}} \approx \sigma_d^{\text{even}}$, where

$$ \sigma_d^{\text{even}} = \sqrt{L^{-1}\partial_{\mu}\langle\varphi\rangle_d^{\text{even}}} = L^{-1/2}\sqrt{\langle\varphi\rangle - 2\langle\varphi\rangle^2}. \qquad (63) $$

The replacement of $D_{N_d}^L$ by $2C_{N_d}^{L/2}$ is exact for the one particle-hole pair state ($N_d = 1$) and at half-filling, ($N_d = L/2$), two cases for which there are no allowed mixed-states (in the latter case because of the high filling of pairs). Although we over-count the (unique) empty state by a factor of 2 (in order to put $\Xi_d^{\text{even}}$ in the free particle form) and under-count the number of particle-hole pair states for $1 < N_d < L/2$, for which there are allowed mixed-states, the physics captured by the approximation $\sigma_d^{\text{even}}$ is sufficiently faithful to that of the original system to be able to account well for the main quantity of interest, namely the standard deviation in the strong anti-cooperativity limit for $L$ even (see Fig. 9).

Inspection of Fig. 3b and Fig. 9 shows that the approximate form $\sigma_d^{\text{even}}$ accounts well for the exact behavior, going to 0 at zero and HF and reaching a maximum at 1/4 filling, not far from the exact maximum position. Because of the approximations inherent in obtaining $\sigma_d^{\text{even}}$, however, the approximate result cannot capture the weak lack of symmetry about 1/4 filling, i.e., the observed skewness towards lower fillings seen in Figs. 3b and 9. This skewness can be seen even more clearly in the PDFs for $L = 14$ and $\langle\varphi\rangle = 1/3$ presented in Fig. 5. It arises because for $1 < N_d < L/2$ the approximation $2C_{N_d}^{L/2}$ underestimates the exact particle-hole pair
<figure><img src="image_11.png" /><figcaption>Figure 9: Comparison between the approximation for the very strong anticooperativity regime $\sigma_d^{\text{even}}$ and the exact result given by the parametric plot with $\langle\varphi\rangle$ and $\sigma$ (see even case in Fig. 3b) for $J = -10$.</figcaption></figure>

multiplicity $D_{N_d}^L$ more severely for lower fillings than for higher ones (below HF). The reason being that at higher fillings it becomes more and more difficult to have allowed mixed states.

The roughly Gaussian shape for the PDF with $L = 14$ presented in Fig. 5 for $\langle\varphi\rangle = 1/3$ can be understood using the approximate particle-hole pair approach that maps the problem onto free quasi-particles (particle-hole pairs). In this case we can adapt the results for the $J = 0$ case to estimate the particle-hole pair PDF for $J \to -\infty$ and $L$ even,

$$P_d^{\text{even}}(N; \langle\varphi\rangle) \equiv \frac{C_{N_d}^{L/2} e^{\mu_d^{\text{even}} N_d}}{\Xi_d^{\text{even}}}, \qquad (64)$$

using the Gaussian approximation:

$$P_d^{\text{even}}(N; \langle\varphi\rangle) \approx \frac{\exp\left[-\frac{(\varphi - \langle\varphi\rangle)^2}{2(\sigma_d^{\text{even}})^2}\right]}{(2\pi)^{1/2}\sigma_d^{\text{even}}L}, \qquad (65)$$

where we have used

$$\Xi_d^{\text{even}}(\langle\varphi\rangle) = \frac{2}{(1 - 2\langle\varphi\rangle)^{L/2}}, \qquad (66)$$

and

$$\mu_d^{\text{even}}(\langle\varphi\rangle) = \ln\left(\frac{2\langle\varphi\rangle}{1 - 2\langle\varphi\rangle}\right). \qquad (67)$$

The Gaussian PDF approximation is very accurate compared with $P_d^{\text{even}}(N; \langle\varphi\rangle)$ for $L = 14$ and $\langle\varphi\rangle = 1/3$, although $P_d^{\text{even}}(N; \langle\varphi\rangle)$ itself underestimates the exact maximum by about 30% and cannot capture the skewness observed in the exact results (see Fig. 5). The particle-hole pair approximation does reproduce accurately, however, the width of the exact PDF.

In the strong anti-cooperativity limit ($J \to -\infty$) for $L$ odd, the system cannot be approximated by non-overlapping particle-hole pairs at HF. We need to account for a defect in order to approach HF. The restricted grand partition function for $0 \le \langle\varphi\rangle \le 1/2$ in this particle-hole pair approximation is

$$\Xi_{-\infty}^{\text{odd}} = \sum_{N_d=0}^{(L-3)/2} O_{N_d}^L e^{\mu N_d} + L e^{\mu(L-1)/2} + L e^{-|J|+\mu(L+1)/2}, \qquad (68)$$
where $(L-1)/2$ is the maximum number of non-overlapping pairs on a lattice of odd size $L$ and the multiplicity $O_{N_d}^L$ counts the number of distinct ways of putting $N_d$ pairs on a periodic lattice of odd size $L$. The last two terms are the defect contributions, an extra hole or particle not making up a particle-hole pair that costs an energy $-J$ for an extra particle and that leads to a multiplicity $L$ counting the number of ways to place an extra hole or particle on a lattice made up of $(L-1)/2$ pairs in sequence.

Since calculating $O_{N_d}^L$ is a non-trivial combinatorial problem, we follow a procedure similar to the one used in the even case and simplify the problem by considering two disjoint lattices of size $(L+1)/2$. If we neglect mixed particle-hole pair microstates for $0 \le N_d \le (L-3)/2$, i.e., those that have pairs on both lattices, we obtain the following approximate restricted particle-hole pair grand partition function

$$ \Xi_d^{\text{odd}} = \sum_{N_d=0}^{(L-3)/2} 2C_{N_d}^{(L+1)/2}e^{\mu N_d} + Le^{\mu(L-1)/2} + Le^{-|J|+\mu(L+1)/2}, \quad (69) $$

which can be evaluated by completing the binomial sum and then subtracting the added terms:

$$ \Xi_d^{\text{odd}} = 2(1+e^{\mu})^{(L+1)/2} - e^{\mu(L-1)/2} - e^{\mu(L+1)/2}(2-Le^{-|J|}) \quad (70) $$

where $C_{N_d}^{(L+1)/2}$ counts the number of distinct ways of placing $N_d$ particles on a lattice of size $(L+1)/2$. $\Xi_d^{\text{odd}}$ correctly counts the multiplicity of the $(L \pm 1)/2$ states by construction. Although this approximation overcounts the one particle-hole pair state by 1 ($L+1$ instead of $L$) and overcounts the empty state by a factor of 2, we believe that it captures the essential physics of the strong anti-cooperativity case for $L$ odd. The above approximation, $\Xi_d^{\text{odd}}$, can be used to calculate $\langle\varphi\rangle_d^{\text{odd}} = (L\Xi_d^{\text{odd}})^{-1}\partial_\mu\Xi_d^{\text{odd}}$ and $\sigma_d^{\text{odd}} = \sqrt{L^{-1}\partial_\mu\langle\varphi\rangle_d^{\text{odd}}}$, an approximate results that could then be compared with the exact ones. We will not, however, pursue this approach any further here.

## 4.3 Infinite range (IR) model

The mapping between the SR and IR Lattice Gas can be established by performing a perturbative cumulant expansion of the grand partition function:

$$ \begin{aligned} \Xi &= \sum_{\{\varphi_i\}} e^{-H_0 - H_{\text{int}}} \\ &= \Xi_0 \langle e^{-H_{\text{int}}} \rangle_0 \\ &\approx \Xi_0 e^{-\langle H_{\text{int}} \rangle_0}, \end{aligned} \quad (71) $$

where the last line gives the first order cumulant expansion. The subscript 0 indicates a statistical average with respect to the non-interacting (Hill-Langmuir) model (ideal Lattice Gas). By calculating $\langle H_{\text{int}} \rangle_0$, we obtain $\langle H_{\text{int}} \rangle_0 = -J_{\text{SR}}L\langle\varphi\rangle_0^2$ for the SR model and $-\frac{1}{2}J_{\text{IR}}L(L-1)\langle\varphi\rangle_0^2$ for the IR model, results that suggest the following mapping between the two models:

$$ J_{\text{SR}} \leftrightarrow \frac{1}{2} J_{\text{IR}}(L-1). \quad (72) $$

This mapping is exact to lowest order in the coupling constant.

We adopt the following strategy to avoid carrying out an explicit calculation of the standard deviation $\sigma$ directly within the cumulant approximation, which is cumbersome. The weak coupling results for the SR model can be found by expanding the exact parametric results for $\langle\varphi\rangle(\mu)$ and the variance $v(\mu) = \sigma^2(\mu)$ in powers of $J$. To first order in $J_{\text{SR}}$ we find

$$ \left(\frac{\sigma_0}{\sigma_{\text{SR}}}\right)^2 = 1 - 2J_{\text{SR}}L\sigma_0^2 + \mathcal{O}(J_{\text{SR}}^2) = 1 - 2J_{\text{SR}}\sigma_\infty^2 + \mathcal{O}(J_{\text{SR}}^2). \quad (73) $$
We can then use the above mapping between the SR and LR models to get a result that is exact to first order in $J_{IR}$ for the IR model:

$$ \left( \frac{\sigma_0}{\sigma_{IR}} \right)^2 = 1 - L(L-1)J_{IR}\sigma_0^2 + \mathcal{O}(J_{IR}^2) = 1 - (L-1)J_{IR}\sigma_\infty^2 + \mathcal{O}(J_{IR}^2). \quad (74) $$

One can check the above method for SR model at half-filling:

$$ \left[ 1 - \left( \frac{\sigma_0}{\sigma_{SR}} \right)^2 \right] = \frac{1}{2} J_{SR} + \mathcal{O}(J_{SR}^2) \quad (\text{half-filling}). \quad (75) $$

It would be useful to find a good interpolation formula for $\sigma$ as an explicit function of $L$, $J$, and $\langle\varphi\rangle$ for both the SR and IR models. A reasonably good interpolation scheme can be set up by recognizing that the exact result for the standard deviation at half-filling for small $J$ is to a very good approximation,

$$ \left( \frac{\sigma_0}{\sigma_{SR}} \right)^2 = \exp(-2J_{SR}\sigma_\infty^2) \times \left[ 1 + \mathcal{O}\left( \left[ \sigma_\infty^2 J_{SR} \right]^L \right) \right] \quad (\text{half-filling}) \quad (76) $$

(where $\sigma_\infty^2 = 1/4$ at HF), which shows that it's best to use an exponential resummation of the small $J$ perturbation expansion.

Comparing the exact SR results with the above approximation shows that an even better approximation is to take the 5/2 root mean (a heuristic choice),

$$ \left( \frac{\sigma_0}{\sigma_{SR}} \right)^2 \approx \left[ \exp(-5J_{SR}\sigma_\infty^2) + L^{-5/2} \right]^{2/5}, \quad (77) $$

to get a smooth interpolation. The above 5/2 root mean square average leads to our best approximation for $\sigma_{SR}$:

$$ \sigma_{SR} \approx \frac{\sigma_0}{[\exp(-5J_{SR}\sigma_\infty^2) + L^{-5/2}]^{1/5}}. \quad (78) $$

This is a convenient form because the dependence on $\langle\varphi\rangle$ (through $\sigma_\infty$), $J_{SR}$, and $L$ is made explicit (since $\sigma_0 = L^{-1/2}\sigma_\infty$). The crossover to the large $J$ saturation limit is centered on

$$ J_{co} \equiv \ln(L)/(2\sigma_{\infty}^{2}). \quad (79) $$

$J_{co}$ depends on $L$ and filling $\langle\varphi\rangle$, increasing slowly with $L$ and reaching a minimum at half-filling for fixed $L$. The convergence of $\sigma$ as a function of $\langle\varphi\rangle$ to the large $J$ saturation limit is uniform only for $J < J_{co}$.

The approximation (78) is very accurate for all fillings (except possibly very close to zero and full-filling) as long as $L^{-5/2} \ll 1$ (the case for $L \approx 10$) and $J < J_{co}$ (see Fig. 16). For $L^{-5/2} \ll 1$ and $J > J_{co}$, this approximation is accurate only within a window centered on half-filling. Because of the factor $\sigma_\infty^2$ in the argument of the exponential, convergence to $\sigma_\infty$ is fastest at HF and becomes slower and slower as moves away from HF to zero and full filling. For $L \approx 10$, the window of accuracy for $J < 2$ is $0.2 < \langle\varphi\rangle < 0.8$, which encompasses the experimental data window.

The interpretation of the above results is in accordance with our previous understanding: as long as the system is far enough from zero and full-filling, increasing $J$ leads to a flattening of the Gaussian PDF with $\sigma$ growing exponentially with $J$ until the system reaches the crossover region, $J \approx J_{co}$, before saturating at $\sigma_\infty$ for $J \gg J_{co}$. When the system reaches the crossover region the Gaussian form is no longer accurate because one (or two) of the wings starts to touch one (or two) of the boundaries, before morphing into a bimodal (two-state) system.

The corresponding result for the IR case can be obtained using the mapping (72). We note in passing that an approximate Gaussian effective free energy approach can be developed for
the IR model along the lines of the one developed above for the non-interacting (Hill-Langmuir) model. This result can be used to find an extended Gaussian approximation for the IR model,

$$P_I R(N; \langle \varphi \rangle) \approx \frac{\exp \left[ - \frac{(\varphi - \langle \varphi \rangle)^2}{2\sigma_{\text{IR}}^2} \right]}{(2\pi)^{1/2}\sigma_{\text{IR}}L}, \quad (80)$$

where the variance $\sigma_{\text{IR}}$ can be approximated using the approximation for the SR model (78) and the mapping (72).

## 4.4 System size, fluctuations, and experimental precision

By proposing a pragmatic operational definition, we clarify here what we mean by a small system. Before doing so we note that if we were to attempt to qualify the size of an adsorption system by the ratio between the number of binding sites, $L$, and the size of a correlated domain (twice the correlation length), we would only be characterizing the Hill-Langmuir nature of the fluctuations, occurring when $2\xi/L \ll 1$, or the bimodal nature of the fluctuations, occurring when $2\xi/L \gg 1$, and not the feasibility of extracting the cooperativity from the experimentally measured relative occupancy fluctuations at equilibrium.

We propose here that a system be considered a small cooperative one if fluctuations, as quantified by the standard deviation of relative occupancy, are within experimental resolution and can be used to extract, using the theoretical framework that we have developed, the amplitude of the interaction strength, $J$, and that this strength be within the usual biophysical range, $0 \le J \le J_{\text{max}} \approx 10$, for adsorption systems.

To quantify our capacity to resolve $J$ theoretically, we therefore propose the following criterion: the derivative of $\sigma$ with respect to $J$, $\sigma'_{\text{HF}}(J; L)$, must be larger than a given value $\Delta\sigma$, whose choice depends on the measurement precision, over a $J$ range starting at 0. In what follows, for the application to the BFM, we shall choose $\Delta\sigma = 0.02$, which is compatible with the experimental error for the data provided in Section 2.5. This choice is motivated by comparing the parametric plots in Fig. 3 to the experimental data (see Fig. 6). One can estimate the coupling $J$ from the data when these fall in regions of the parametric plots where the lines corresponding to different values of $J$ are spaced far enough apart compared to the experimental error. The highest resolution for positive values of $J$, and fixed $L \ge 3$, is achieved at half-filling (HF), and for a value $J^* \approx 2\ln(L/2)$ that maximizes the derivative of the standard deviation. These maxima correspond to the peaks of the curves in Fig. 10a where we plot $\sigma'_{\text{HF}}(J; L)$ as a function of $J$ for different but fixed values of $L$ and compare it with the threshold value $\Delta\sigma = 0.02$.

We observe that for the choice of $\Delta\sigma = 0.02$, $L^* \approx 39$ is the maximum system size for which there is a feasibility window ($0 < J < 8$) starting at 0 for extracting the cooperativity from the experimental data. This threshold value $L^*$ can be recovered analytically by using the approximation Eq. (78) for the standard deviation, which is highly accurate near half-filling provided $L \gg 1$, to find $\sigma'_{\text{HF}}(J = 0; L) \approx 1/(8L^{1/2})$. One can then solve $\sigma'_{\text{HF}}(J = 0; L) = \Delta\sigma = 0.02$ to obtain $L^* = 1/(64\Delta\sigma^2)$. The same approximation leads to $\sigma''_{\text{HF}}(J = 0; L) \approx 1/(32L^{1/2})$. These approximations show that both the first and second derivatives of $\sigma_{\text{HF}}$ at zero coupling ($J = 0$) decrease as $L^{-1/2}$. For high coupling, $J \gg J_{\text{max}}$, the system assumes a bimodal character and the standard deviation saturates at the infinite coupling limit.

We conclude, based on our definition and our choice of $\Delta\sigma$, that a system is small if $L \le 39$. We also observe that for any system of finite size with $L > 39$, there is always a feasibility window for sufficiently strongly cooperative systems, $J_{\text{min}} < J < J_{\text{max}}$, characterized by a peak at a value of $J = J^*$ that increases logarithmically with increasing $L$, but with $J_{\text{min}} > 0$. Consequently, for $L \gg 39$, the standard deviation at HF is insensitive to $J$ both at low coupling, where the system acts as a Hill-Langmuir one, and high coupling, where the system becomes bimodal.
<figure><img src="image_12.png" /><figcaption>Figure 10: (a) Derivative of the standard deviation at half-filling, $\sigma'_{HF}(J; L) \equiv d\sigma_{HF}/dJ$, plotted for systems of size $L = 10, 39, 10^2, 10^3, 10^4$, respectively. The horizontal dashed gray line represents the threshold value $\Delta\sigma = 0.02$. (b) Contour plot of $\sigma'_{HF}(J; L)$ for positive values of $J$ and $L \ge 3$. The red contour lines correspond to $\sigma'_{HF}(J; L) = \Delta\sigma = 0.02$.</figcaption></figure>

The (semi-log) contour plot in Fig. 10b gives the isolines of $\sigma'_{HF}(J)$ as a function of the interaction potential $J$ and the system size $L$, where the latter, for the sake of simplicity, is treated as a continuous variable. Recall the asymptotic behavior of the correlation length at half-filling for large values of the interaction potential [19]: $J \sim 2 \ln(2\xi_{HF})$. Hence, if we were to plot the logarithm of $\xi_{HF}$ instead of $J$ on the horizontal axis, we would get a qualitatively equivalent contour diagram. The two red contour lines correspond to the constant value $\Delta\sigma = 0.02$. The region between them indicates the set of parameters $J$ and $L$ for which $\sigma'_{HF}(J) > \Delta\sigma$. We finally conclude that a system of size $L = 13$, like the one considered in the main text to describe the BFM, is well within the feasibility window for small systems.

## 4.5 Effective Hill coefficient

Following the analysis in [7] on the logarithmic sensitivity of kinetics schemes, we can also quantify cooperative behavior by introducing an effective Hill coefficient which, for our system, takes the following form:

$$H_{\text{eff}} = \frac{d}{d\mu} \ln \left( \frac{\langle \varphi \rangle}{1 - \langle \varphi \rangle} \right) = \frac{\sigma^2 L}{\langle \varphi \rangle (1 - \langle \varphi \rangle)}. \quad (81)$$

Here, we have exploited the relation $\partial_{\mu}\langle\varphi\rangle = \sigma^2 L$.

The effective Hill function defined above takes on a maximal value $H_{\text{eff}}^*$ at half-filling:

$$H_{\text{eff}}^*(J) = 4\sigma_{\text{HF}}^2 L, \quad (82)$$

where $\sigma_{\text{HF}}^2$ is the square of the relative standard deviation at half-filling whose explicit dependency on $J$ is captured by Eq. (5). This expression tends rapidly to zero for negative values of $J$ (anticooperativity). Moreover, for $J = 0$, the coefficient equals 1, as expected for the Hill-Langmuir case. Finally, by increasing the interaction potential, $\sigma_{\text{HF}}$ tends to $1/2$, and so $H_{\text{eff}}^* \to L$. This result is in agreement with the rule that at equilibrium, the Hill coefficient is expected not to exceed the maximal number of ligands that can bind simultaneously, which, in the case of an adsorption process, corresponds to the number of binding sites, $L$. Fig. 11 shows the sigmoidal allure of the Hill coefficient, bounded from above by the system size and below by 0.

It should be noted, however, that in [7] it has been demonstrated that this upper bound for the effective Hill coefficient may be surpassed when the system is out-of-equilibrium. In such
<figure><img src="image_13.png" /><figcaption>Figure 11: The effective Hill coefficient at half-filling from Eq. (82) as a function of the interaction potential *J* for a system of size *L* = 13.</figcaption></figure>

cases, the new limit is determined not by the size of the lattice, but by the number of states of
the kinetic scheme whose exit rates increase with the concentration of the ligands.

4.6 Supplemental plots

<figure><img src="image_14.png" /><figcaption>Figure 12: Contour plot of the mean relative occupancy in equilibrium, $\langle\phi\rangle$, as a function of the dimensionless interaction potential $J$ and chemical potential $\mu$ for a system of size $L = 14$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13: Contour plot of the standard deviation of the relative occupancy at equilibrium, $\sigma = \sqrt{\langle\varphi^2\rangle - \langle\varphi\rangle^2}$, as a function of the dimensionless interaction potential $J$ and chemical potential $\mu$ for a system of size $L = 14$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14: Standard deviation of the mean occupancy at half-filling, $\langle\varphi\rangle = 1/2$, as a function of the dimensionless interaction potential $J$ for a system of size $L = 13$ (blue) and $L = 14$ (orange). For $J \ge 0$ the prediction of the block domain approximation, Eq. (13), is indistinguishable from the exact result.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Correlation length, $\xi$, in the limit of $J \to \infty$. It diverges at half-filling and is mainly higher than $L/2$ (indicated by the red line), except near zero and full fillings.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16: Comparison between the exact result for the standard deviation (solid lines) and the approximate expression (dashed lines) that depends explicitly on $\langle\phi\rangle$, $J$, and $L$ (Eq. (78)).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17: Fit of the theoretical standard deviation (blue solid line) as a function of $\langle\phi\rangle$ (for $L=13$) to the BFM experimental data (blue data points with experimental error bars). The green shaded area is the 90% confidence window and the best fit value is $J = 1.21 \pm 0.22$, where 0.22 is the standard error of the nonlinear fit. The solid purple curve is the result for the Hill-Langmuir model ($J=0$) and the solid black curve is the strong coupling result ($J \to \infty$).</figcaption></figure>
## 4.7 Experimental data

<table><caption>Table 1: Data of average occupancy and standard deviation corresponding to Fig. 6. Data are obtained by averaging over several traces [40, 55]. As there are not enough traces for good statistics (from 30 to 40), data is also averaged over several time points at equilibrium.</caption><thead><tr><th>Bead size (drag)</th><th>⟨φ⟩</th><th>σ</th></tr></thead><tbody><tr><td rowspan="3">300 nm</td><td>0.29 ± 0.03</td><td>0.17 ± 0.02</td></tr><tr><td>0.355 ± 0.006</td><td>0.153 ± 0.011</td></tr><tr><td>0.41 ± 0.02</td><td>0.15 ± 0.02</td></tr><tr><td>300 nm + Glycerol</td><td>0.461 ± 0.005</td><td>0.175 ± 0.007</td></tr><tr><td rowspan="3">500 nm</td><td>0.56 ± 0.02</td><td>0.18 ± 0.02</td></tr><tr><td>0.592 ± 0.011</td><td>0.21 ± 0.02</td></tr><tr><td>0.67 ± 0.02</td><td>0.18 ± 0.04</td></tr><tr><td>500 nm + Glycerol</td><td>0.705 ± 0.005</td><td>0.158 ± 0.005</td></tr><tr><td rowspan="4">1300 nm</td><td>0.74 ± 0.03</td><td>0.13 ± 0.02</td></tr><tr><td>0.80 ± 0.03</td><td>0.153 ± 0.012</td></tr><tr><td>0.828 ± 0.006</td><td>0.18 ± 0.04</td></tr><tr><td>0.826 ± 0.011</td><td>0.161 ± 0.009</td></tr></tbody></table>

<table><caption>Table 2: Data of probability distribution corresponding to Fig.7. Data are obtained by averaging over several time points at equilibrium the state of all traces.</caption><thead><tr><th rowspan="3"># Stator units</th><th colspan="3">Bead size</th></tr><tr><th>300 nm</th><th>500 nm</th><th>1300 nm</th></tr><tr><th colspan="3">Probability</th></tr></thead><tbody><tr><td>0</td><td>0.12 ± 0.04</td><td>0.005 ± 0.014</td><td>0.01 ± 0.20</td></tr><tr><td>1</td><td>0.09 ± 0.04</td><td>0.005 ± 0.013</td><td>0</td></tr><tr><td>2</td><td>0.11 ± 0.05</td><td>0</td><td>0</td></tr><tr><td>3</td><td>0.15 ± 0.04</td><td>0.01 ± 0.02</td><td>0</td></tr><tr><td>4</td><td>0.15 ± 0.04</td><td>0.05 ± 0.02</td><td>0</td></tr><tr><td>5</td><td>0.22 ± 0.04</td><td>0.04 ± 0.03</td><td>0</td></tr><tr><td>6</td><td>0.13 ± 0.03</td><td>0.17 ± 0.05</td><td>0</td></tr><tr><td>7</td><td>0.033 ± 0.014</td><td>0.29 ± 0.10</td><td>0.002 ± 0.011</td></tr><tr><td>8</td><td>0</td><td>0.13 ± 0.04</td><td>0.12 ± 0.04</td></tr><tr><td>9</td><td>0</td><td>0.13 ± 0.04</td><td>0.14 ± 0.03</td></tr><tr><td>10</td><td>0</td><td>0.05 ± 0.02</td><td>0.28 ± 0.07</td></tr><tr><td>11</td><td>0</td><td>0.04 ± 0.02</td><td>0.25 ± 0.06</td></tr><tr><td>12</td><td>0</td><td>0.06 ± 0.02</td><td>0.18 ± 0.05</td></tr><tr><td>13</td><td>0</td><td>0.03 ± 0.02</td><td>0.02 ± 0.03</td></tr></tbody></table>
<table><caption>Table 3</caption><thead><tr><th rowspan="2">Bead size</th><th colspan="3">μ estimated from &lang;φ&rang;</th></tr><tr><th>J = 1</th><th>J = 2</th><th>J = 3</th></tr></thead><tbody><tr><td rowspan="3">300 nm</td><td>-1.56 ± 0.08</td><td>-2.35 ± 0.05</td><td>-3.21 ± 0.03</td></tr><tr><td>-1.37 ± 0.02</td><td>-2.222 ± 0.011</td><td>-3.14 ± 0.01</td></tr><tr><td>-1.23 ± 0.04</td><td>-2.14 ± 0.03</td><td>-3.08 ± 0.02</td></tr><tr><td>300 nm + Glycerol</td><td>-1.096 ± 0.012</td><td>-2.058 ± 0.007</td><td>-3.035 ± 0.004</td></tr><tr><td rowspan="3">500 nm</td><td>-0.85 ± 0.05</td><td>-1.91 ± 0.03</td><td>-2.95 ± 0.02</td></tr><tr><td>-0.77 ± 0.03</td><td>-1.86 ± 0.02</td><td>-2.92 ± 0.01</td></tr><tr><td>-0.56 ± 0.06</td><td>-1.73 ± 0.04</td><td>-2.84 ± 0.02</td></tr><tr><td>500 nm + Glycerol</td><td>-0.46 ± 0.02</td><td>-1.67 ± 0.01</td><td>-2.80 ± 0.01</td></tr><tr><td rowspan="3">1300 nm</td><td>-0.32 ± 0.11</td><td>-1.58 ± 0.07</td><td>-2.75 ± 0.04</td></tr><tr><td>-0.12 ± 0.011</td><td>-1.45 ± 0.07</td><td>-2.67 ± 0.05</td></tr><tr><td>0.012 ± 0.032</td><td>-1.37 ± 0.02</td><td>-2.614 ± 0.013</td></tr></tbody></table>

<table><caption>Table 4</caption><thead><tr><th rowspan="2">Bead size</th><th colspan="3">μ estimated from σ</th></tr><tr><th>J = 1</th><th>J = 2</th><th>J = 3</th></tr></thead><tbody><tr><td rowspan="3">300 nm</td><td>-1.4 ± 0.4</td><td>-2.55 ± 0.11</td><td>-3.48 ± 0.06</td></tr><tr><td>-1.61 ± 0.21</td><td>-2.63 ± 0.09</td><td>-3.53 ± 0.05</td></tr><tr><td>-1.59 ± 0.33</td><td>-2.6 ± 1.1</td><td>-3.53 ± 0.07</td></tr><tr><td>300 nm + Glycerol</td><td>-1.19 ± 0.19</td><td>-1.5 ± 0.9</td><td>-3.45 ± 0.3</td></tr><tr><td rowspan="3">500 nm</td><td>-0.99 ± 0.01</td><td>-1.53 ± 0.13</td><td>-3.44 ± 0.06</td></tr><tr><td>-1.00 ± 0.01</td><td>-2.30 ± 0.12</td><td>-3.4 ± 0.7</td></tr><tr><td>-0.99 ± 0.01</td><td>-1.6 ± 0.5</td><td>-3.4 ± 0.7</td></tr><tr><td>500 nm + Glycerol</td><td>-1.53 ± 0.09</td><td>-2.60 ± 0.03</td><td>-3.51 ± 0.02</td></tr><tr><td rowspan="3">1300 nm</td><td>-0.2 ± 1.4</td><td>-2.75 ± 0.16</td><td>-3.6 ± 0.1</td></tr><tr><td>-1.61 ± 0.20</td><td>-2.63 ± 0.08</td><td>-3.53 ± 0.05</td></tr><tr><td>-1.00 ± 0.01</td><td>-1.6 ± 0.5</td><td>-3.4 ± 0.7</td></tr></tbody></table>
References

[1] A. Whitty, "Cooperativity and biological complexity," *Nature chemical biology*, vol. 4, no. 8, pp. 435-439, 2008.

[2] R. T. Ogata and H. M. McConnell, "Mechanism of cooperative oxygen binding to hemoglobin," *Proceedings of the National Academy of Sciences*, vol. 69, no. 2, pp. 335-339, 1972.

[3] K. A. Dill, K. M. Fiebig, and H. S. Chan, "Cooperativity in protein-folding kinetics.," *Proceedings of the National Academy of Sciences*, vol. 90, no. 5, pp. 1942-1946, 1993.

[4] K. J. Tielrooij, N. Garcia-Araez, M. Bonn, and H. J. Bakker, "Cooperativity in ion hydration," *Science*, vol. 328, no. 5981, pp. 1006-1009, 2010.

[5] J. R. Williamson, "Cooperativity in macromolecular assembly," *Nature chemical biology*, vol. 4, no. 8, pp. 458-465, 2008.

[6] S. A. West, S. P. Diggle, A. Buckling, A. Gardner, and A. S. Griffin, "The social lives of microbes," *Annual Review of Ecology, Evolution, and Systematics*, vol. 38, pp. 53-77, 2007.

[7] J. A. Owen and J. M. Horowitz, "Size limits the sensitivity of kinetic schemes," *Nature Communications*, vol. 14, no. 1, p. 1280, 2023.

[8] C. H. Berg, "The rotary motor of bacterial flagella," *Annual Review of Biochemistry*, vol. 72, no. 1, pp. 19-54, 2003.

[9] J. A. Nirody, Y.-R. Sun, and C.-J. Lo, "The biophysicist's guide to the bacterial flagellar motor," *Advances in Physics: X*, vol. 2, no. 2, pp. 324-343, 2017.

[10] K. I. Ito, S. Nakamura, and S. Toyabe, "Cooperative stator assembly of bacterial flagellar motor mediated by rotation," *Nature Communications*, vol. 12, no. 3218, 2021.

[11] F. Bai, R. W. Branch, D. V. Nicolau, T. Pilizota, B. C. Steel, P. K. Maini, and R. M. Berry, "Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch," *Science*, vol. 327, no. 5966, pp. 685-689, 2010.

[12] H. Fukuoka, T. Wada, S. Kojima, A. Ishijima, and M. Homma, "Sodium-dependent dynamic assembly of membrane complexes in sodium-driven flagellar motors," *Molecular microbiology*, vol. 71, pp. 825-35, 02 2009.

[13] M. J. Tipping, N. J. Delalez, R. Lim, R. M. Berry, and J. P. Armitage, "Load-dependent assembly of the bacterial flagellar motor.," *mBio*, vol. 4, no. 4, 2013.

[14] D. Castillo, S. Nakamura, Y. Morimoto, Y.-S. Che, N. Kami-ike, S. Kudo, T. Minamino, and K. Namba, "The c-terminal periplasmic domain of motb is responsible for load-dependent control of the number of stators of the bacterial flagellar motor," *BIOPHYSICS*, vol. 9, pp. 173-181, 12 2013.

[15] N. Wadhwa, R. Phillips, and H. C. Berg, "Torque-dependent remodeling of the bacterial flagellar motor.," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 116, no. 24, pp. 11764-11769, 2019.

[16] J. A. Nirody, A. L. Nord, and R. M. Berry, "Load-dependent adaptation near zero load in the bacterial flagellar motor.," *Journal of the Royal Society, Interface*, vol. 16, no. 159, p. 20190300, 2019.
[17] F. Bai, T. Minamino, Z. Wu, K. Namba, and J. Xing, "Coupling between switching regulation and torque generation in bacterial flagellar motor," *Phys. Rev. Lett.*, vol. 108, p. 178105, Apr 2012.

[18] B. M. McCoy and T. T. Wu, *The two-dimensional Ising model*. Courier Corporation, 2014.

[19] C. Godreche and J. Luck, "Response of non-equilibrium systems at criticality: exact results for the glauber-ising chain," *Journal of Physics A: Mathematical and General*, vol. 33, no. 6, p. 1151, 2000.

[20] E. M. Purcell, "Life at low reynolds number," *American Journal of Physics*, vol. 45, no. 1, pp. 3-11, 1977.

[21] D. Kearns, "A field guide to bacterial swarming motility," *Nat Rev Microbiol*, vol. 8, p. 634-644, 2010.

[22] A. M. Spormann, "Gliding motility in bacteria: Insights from studies of myxococcus xanthus," *Microbiology and Molecular Biology Reviews*, vol. 63, no. 3, pp. 621-641, 1999.

[23] J. Henrichsen, "Twitching motility," *Annual Review of Microbiology*, vol. 37, no. 1, pp. 81-93, 1983.

[24] A. Martínez, S. Torello, and R. Kolter, "Sliding motility in mycobacteria," *Journal of Bacteriology*, vol. 181, no. 23, pp. 7331-7338, 1999.

[25] M. Dvorashyna and E. Lauga, "Hydrodynamics and direction change of tumbling bacteria," *PLOS ONE*, vol. 16, pp. 1-26, 07 2021.

[26] N. Wadhwa and H. Berg, "Bacterial motility: machinery and mechanisms," *Nat Rev Microbiol*, vol. 20, p. 161-173, 2022.

[27] S. Nakamura and T. Minamino, "Flagella-driven motility of bacteria," *Biomolecules*, vol. 9, no. 7, 2019.

[28] M. Grognot and K. M. Taute, "More than propellers: how flagella shape bacterial motility behaviors," *Current Opinion in Microbiology*, vol. 61, pp. 73-81, 2021.

[29] M. T. Brown, B. C. Steel, C. Silvestrin, D. A. Wilkinson, N. J. Delalez, C. N. Lumb, B. Obara, J. P. Armitage, and R. M. Berry, "Flagellar hook flexibility is essential for bundle formation in swimming escherichia coli cells.," *Journal of bacteriology*, vol. 194, no. 13, pp. 3495-3501, 2012.

[30] C. V. Gabel and H. C. Berg, "Speed of the flagellar rotary motor of Escherichia coli varies linearly with protonmotive force," *Proc Natl Acad Sci U S A*, vol. 100, no. 15, pp. 8748-51, 2003.

[31] E. R. Hosking, C. Vogt, E. P. Bakker, and M. D. Manson, "Escherichia coli MotAB proton channel unplugged," *Journal of Molecular Biology*, vol. 364, no. 5, pp. 921-937, 2006.

[32] Y. Sowa, M. Homma, A. Ishijima, and R. M. Berry, "Hybrid-fuel bacterial flagellar motors in escherichia coli.," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 111, no. 9, pp. 3436-3441, 2014.

[33] W. S. Ryu, R. M. Berry, and H. C. Berg, "Torque-generating units of the flagellar motor of escherichia coli have a high duty ratio," *Nature*, vol. 403, no. 6768, pp. 444-7, 2000.
[34] S. W. Reid, M. C. Leake, J. H. Chandler, C. J. Lo, J. P. Armitage, and R. M. Berry, "The maximum number of torque-generating units in the flagellar motor of escherichia coli is at least 11," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 103, no. 21, pp. 8066-71, 2006.

[35] D. C. Fung and H. C. Berg, "Powering the flagellar motor of escherichia coli with an external voltage source," *Nature*, vol. 375, pp. 809-812, 1995.

[36] S. Kojima, N. Nonoyama, N. Takekawa, H. Fukuoka, and M. Homma, "Mutations targeting the C-terminal domain of FliG can disrupt motor assembly in the Na(+)-driven flagella of Vibrio alginolyticus.," *Journal of Molecular Biology*, vol. 414, no. 1, pp. 62-74, 2011.

[37] A. Biquet-Bisquert, G. Labesse, F. Pedaci, and A. L. Nord, "The dynamic ion motive force powering the bacterial flagellar motor," *Frontiers in Microbiology*, vol. 12, no. 659464, 2021.

[38] P. P. Lele, R. W. Branch, V. S. J. Nathan, and H. C. Berg, "Mechanism for adaptive remodeling of the bacterial flagellar switch," *Proceedings of the National Academy of Sciences*, vol. 109, no. 49, pp. 20018-20022, 2012.

[39] P. P. Lele, B. G. Hosu, and H. C. Berg, "Dynamics of mechanosensing in the bacterial flagellar motor.," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 110, no. 29, pp. 11839-44, 2013.

[40] A. L. Nord, E. Gachon, R. Perez-Carrasco, J. A. Nirody, A. Barducci, R. M. Berry, and F. Pedaci, "Catch bond drives stator mechanosensitivity in the bacterial flagellar motor," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 114, no. 49, pp. 12952-12957, 2017.

[41] X. Hu, F. M. Margadant, M. Yao, and M. P. Sheetz, "Molecular stretching modulates mechanosensing pathways," *Protein Science*, vol. 26, no. 7, pp. 1337-1351, 2017.

[42] N. Wadhwa, Y. Tu, and H. C. Berg, "Mechanosensitive remodeling of the bacterial flagellar motor is independent of direction of rotation," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 118, no. 15, p. e2024608118, 2021.

[43] V. Sourjik and N. S. Wingreen, "Responding to chemical gradients: bacterial chemotaxis," *Current Opinion in Cell Biology*, vol. 24, no. 2, pp. 262-268, 2012. Cell regulation.

[44] Q. Ma, Y. Sowa, M. Baker, and F. Bai, "Bacterial flagellar motor switch in response to chey-p regulation and motor structural alterations," *Biophysical Journal*, vol. 110, pp. 1411-1420, 03 2016.

[45] J. Yuan and H. Berg, "Ultrasensitivity of an adaptive bacterial motor," *Journal of molecular biology*, vol. 425, 02 2013.

[46] S. E. Tusk, N. J. Delalez, and R. M. Berry, "Subunit exchange in protein complexes," *Journal of Molecular Biology*, vol. 430, no. 22, pp. 4557-4579, 2018. Plasticity of Multi-Protein Complexes.

[47] M. Santiveri, A. Roa-Eguiara, C. Kühne, N. Wadhwa, H. Hu, H. C. Berg, M. Erhardt, and N. M. I. Taylor, "Structure and function of stator units of the bacterial flagellar motor," *Cell*, vol. 183, no. 1, pp. 244-257, 2020.

[48] J. Deme, S. Johnson, O. Vickery, A. Muellbauer, H. Monkhouse, T. Griffiths, R. H. James, B. C. Berks, J. W. Coulton, P. Stansfeld, and S. M. Lea, "Structures of the stator complex that drives rotation of the bacterial flagellum," *Nature Microbiology*, vol. 5, p. 1553-1564, 2020.
[49] S. Zhu, M. Takao, N. Li, M. Sakuma, Y. Nishino, M. Homma, S. Kojima, and K. Imada, "Conformational change in the periplamic region of the flagellar stator coupled with the assembly around the rotor.," *Proceedings of the National Academy of Sciences of the United States of America*, vol. 111, no. 37, pp. 13523-13528, 2014.

[50] S. Kojima, K. Imada, M. Sakuma, Y. Sudo, C. Kojima, T. Minamino, M. Homma, and K. Namba, "Stator assembly and activation mechanism of the flagellar motor by the periplasmic region of motb," *Molecular Microbiology*, vol. 73, no. 4, pp. 710-8, 2009.

[51] M. J. Tipping, B. C. Steel, N. J. Delalez, R. M. Berry, and J. P. Armitage, "Quantification of flagellar motor stator dynamics through in vivo proton-motive force control," *Molecular Microbiology*, vol. 87, no. 2, pp. 338-347, 2013.

[52] S. Kojima, M. Takao, G. Almira, I. Kawahara, M. Sakuma, M. Homma, C. Kojima, and K. Imada, "The helix rearrangement in the periplasmic domain of the flagellar stator b subunit activates peptidoglycan binding and ion influx," *Structure*, vol. 26, no. 4, pp. 590-598, 2018.

[53] N. Terahara, N. Kodera, T. Uchihashi, T. Ando, K. Namba, and T. Minamino, "Na-induced structural transition of motps for stator assembly of the flagellar motor.," *Science Advances*, vol. 3, no. 11, p. 4119, 2017.

[54] T.-S. Lin, S. Zhu, S. Kojima, M. Homma, and C.-J. Lo, "Flil association with flagellar stator in the sodium-driven vibrio motor characterized by the fluorescent microscopy.," *Scientific reports*, vol. 8, no. 1, p. 11172, 2018.

[55] R. Perez-Carrasco, M.-J. Franco-Oñate, J.-C. Walter, J. Dorignac, F. Geniet, J. Palmeri, A. Parmeggiani, N.-O. Walliser, and A. L. Nord, "Relaxation time asymmetry in stator dynamics of the bacterial flagellar motor," *Science Advances*, vol. 8, no. 12, p. eabl8112, 2022.

[56] T. Duke, N. Le Novère, and D. Bray, "Conformational spread in a ring of proteins: a stochastic approach to allostery," *Journal of Molecular Biology*, vol. 308, no. 3, pp. 541-553, 2001.

[57] M. Kardar, *Statistical Physics of Fields*. Cambridge University Press, 2007.