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
email: john.palmeri@umontpellier.fr

Nils-Ole Walliser
email: nils-ole.walliser@umontpellier.fr

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

Keywords:

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. Finally, we estimate the stator-stator interaction potential and attempt to quantify the adaptability of the BFM.

Introduction

Cooperativity is a pervasive phenomenon that emerges in biological systems at different time and length scales and levels of complexity [START_REF] Whitty | Cooperativity and biological complexity[END_REF]: from the molecular scale [START_REF] Ogata | Mechanism of cooperative oxygen binding to hemoglobin[END_REF][START_REF] Dill | Cooperativity in protein-folding kinetics[END_REF][START_REF] Tielrooij | Cooperativity in ion hydration[END_REF][START_REF] Williamson | Cooperativity in macromolecular assembly[END_REF] to multiple-cell organization [START_REF] West | The social lives of microbes[END_REF], 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 [START_REF] Owen | Size limits the sensitivity of kinetic schemes[END_REF]. 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 coupled to a thermal heat bath, 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. Finally, we apply the model by comparing our theoretical predictions with experimental data for the bacterial flagellar motor (BFM) of E. coli [START_REF] Berg | The rotary motor of bacterial flagella[END_REF][START_REF] Nirody | The biophysicist's guide to the bacterial flagellar motor[END_REF], a macro-molecular complex that has previously been found to exhibit cooperative binding of the torque-generating (stator) units [START_REF] Ito | Cooperative stator assembly of bacterial flagellar motor mediated by rotation[END_REF]. In this previous work, although evidence for cooperativity was uncovered in the stator unit assembly dynamics, no estimate of the stator-stator interaction potential was proposed, which is our goal here.

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 B T for the short-range model, a result that is coherent with this system's expected smooth adaptive nature to external stimuli [START_REF] Fukuoka | Sodium-dependent dynamic assembly of membrane complexes in sodium-driven flagellar motors[END_REF][START_REF] Tipping | Load-dependent assembly of the bacterial flagellar motor[END_REF][START_REF] Castillo | The c-terminal periplasmic domain of motb is responsible for load-dependent control of the number of stators of the bacterial flagellar motor[END_REF][START_REF] Wadhwa | Torque-dependent remodeling of the bacterial flagellar motor[END_REF][START_REF] Nirody | Load-dependent adaptation near zero load in the bacterial flagellar motor[END_REF].

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.

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 µ. We assign to each binding site i (i = 1, ..., L) a binary variable φ i : φ i = 1 means that a particle (ligand) occupies the i-th site; otherwise, φ i = 0. The array φ = {φ i } = (φ 1 , ..., φ L ) uniquely defines one of the 2 L possible microscopic configurations (or microstates). The relative occupancy in a given microstate is φ = L -1 L i=1 φ 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φ. There are C L N = L!/ [N ! (L -N )!] microstates φ consistent with the occupancy N ; in this case, we say that the mesostate N has multiplicity C L N . The energy of the system depends on the specific microscopic configuration φ according to the grand canonical effective 1D short-range interaction Hamiltonian:

βH(φ) = -J L i=1 φ i φ i+1 -µ L i=1 φ i . (1)
Since we choose periodic boundary conditions, i = L + 1 corresponds to i = 1. Furthermore, we introduce β = 1/(k B T ) with the Boltzmann constant k B , rendering the two adjustable control parameters of our model, J and µ, dimensionless. The effective chemical potential, µ, is considered to be made up of two parts, µ = µ r -ε, where µ r is the chemical potential of the external particle reservoir and ε < 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 µ through µ r and/or ε. When J = 0 we recover the usual noninteracting (Langmuir) model. The predictions for this short-range model will be compared below with those for the longrange model where all pairs of particles interact with the same strength independent of their relative distance on the lattice.

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 [START_REF] Mccoy | The two-dimensional Ising model[END_REF] for a comprehensive overview). It is, therefore, possible to derive analytical expressions for the partition function and 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 partition function Z with respect to µ, we can calculate the relative mean occupancy at equilibrium,

⟨φ⟩ = L -1 L i=1 ⟨φ i ⟩ = (LZ) -1 ∂ µ Z: ⟨φ⟩ = 1 2 1 + tanh L 2ξ sinh X sinh 2 X + e -J , (2)
where X = (J + µ) /2, and ξ is the equilibrium correlation length of the system (given in number of lattice sites), which can be rewritten in terms of J and µ as follows (we refer to Section 4 for more details):

ξ = 1/ ln cosh X + sinh 2 X + e -J cosh X -sinh 2 X + e -J . (3)
The mean equilibrium occupancy equals ⟨N ⟩ = L⟨φ⟩. 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 2) and ( 4), for a system of size L = 13.

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 the relative fluctuations of occupancy are ≳ 1, which encompasses the non-interacting (Langmuir) model for L sufficiently small. Fig. 1a depicts the relative mean occupancy from Eq. ( 2) as a function of the two control parameters of the system: J and µ. The half-filling (HF) contour line ⟨φ⟩ = 0.5 cuts the plot diagonally, meaning that the lattice is, on average, half-filled whenever J = -µ (or X = 0) holds. Below (above) that line, the system's mean occupancy is consistently lower (higher) than 50 %. Contour lines run together close to the plot's bottom right sector diagonal, corresponding to positive J and negative µ. In this sector defined by X ≈ 0, for J ≫ 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 µ. On the contrary, more adaptive systems, that is, systems whose occupancy varies more smoothly from zero to full occupancy as a function of µ, would restrict themselves to the central and upper left sector of Fig. 1a.

For the short-range lattice gas model from Eq. ( 1), the standard deviation σ of the occupancy at equilibrium can be expressed in terms of the derivative with respect to the chemical potential as follows (cf. Section 4): σ = ⟨φ 2 ⟩ -⟨φ⟩ 2 = L -1/2 ∂ µ ⟨φ⟩. Applying this result to Eq. ( 2), we obtain an analytical expression for σ:

σ = 1 2 sech 2 L 2ξ sinh 2 X sinh 2 X + e -J + tanh L 2ξ e -J cosh X L sinh 2 X + e -J 3/2 1/2 , (4)
whose dependency on J and µ is depicted in Fig. 1b. In what follows, we will use σ as a measure of fluctuations, hence the stochastic nature of the system.

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:

σ HF (J) = e J/4 2 √ L tanh L 2ξ HF (J) , (5)
with

ξ HF = 1/ ln 1 + e -J/2 1 -e -J/2 = 1/ ln[coth(J/4)], (6)
which represents the maximum equilibrium correlation length at fixed J. In the Hill-Langmuir limit (J = 0), the correlation length vanishes, and the standard deviation at HF becomes 1/(2 √ L). For low couplings (0 ≤ J < 2), the correlation length increases sub-linearly with J and is accurately given by ξ HF ≈ -1/ ln(J/4). In this same limit, the standard deviation increases exponentially with J, σ HF (J)/σ HF (0) ≈ e J/4 , an approximation valid for J < 4, when L ≈ 10. This increase can be interpreted as a widening of the probability distribution function of states of a given N around the average value ⟨N ⟩, where the probability is a maximum. In the high cooperativity limit J ≫ 1, ξ HF shows an asymptotic exponential growth, ξ HF ≃ 1 2 e J/2 , which is the analog of the well-known zero temperature critical behavior of spin chains [START_REF] Godreche | Response of non-equilibrium systems at criticality: ferromagnetic models in dimension two and above[END_REF]; and consequently, for fixed L, σ HF tends to 1/2 when the correlation length reaches 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, σ HF vanishes as L -1/2 when L → ∞, as expected, in the limit where L is kept fixed, and J tends to +∞ the standard deviation at HF reaches the finite (maximum) value of one-half.

Fluctuations decrease as J < 0 becomes more negative. Along the half-filling diagonal (X = 0) in Fig. 1b, they saturate at a finite value as J → -∞. In this strong anti-cooperative limit, for small enough systems, the result depends sensitively on whether the number of lattice sites is even or odd. For L even, σ 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 (X = 0) can clearly detect the difference between small systems with even and odd numbers of sites.

Fluctuations of the mean occupancy and their dependency on ⟨φ⟩ and J

Often, the experimental parameter one can measure (or control indirectly, as is the case for the BFM, see below) is not the chemical potential of the reservoir µ, but the mean number of bound ligands at equilibrium, ⟨N ⟩. It is, therefore, insightful to plot the standard deviation σ as a function of the mean relative occupancy ⟨φ⟩, which, in general, can be done using a parametric plot.

To unfold this relation, we express ⟨φ⟩ and σ, for different but fixed values of J, as functions of µ and plot this family of curves in Fig. 2a and Fig. 2b for an odd and even number of lattice sites, respectively. Each curve for a given J represents the set of points of coordinates (⟨φ⟩ = f J (µ), σ = g J (µ)) with µ 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: ⟨φ⟩ HF = f J (-J) = 1/2 (as is expected from the underlining particle-hole symmetry of the system under study). The system has zero fluctuations when ⟨φ⟩ = 0 (empty substrate) or ⟨φ⟩ = 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 √ L) for J = 0 and 1/2 for J → +∞, and the σ vs. ⟨φ⟩ curves converge asymptotically towards the limiting J → +∞ one. More subtly, results for the anti-cooperative limit J → -∞ 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 ≈ 10 and mean relative occupancies not too close to 0 or 1, in the range 0 < J < 7 the curves of σ vs. φ 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 ≈ 10 and mean relative occupancies close to 1/2 in the range -3 < J < 0. These results are illustrated in Fig. 2a for L = 13 and Fig. 2b 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 exact or approximate analytical expressions for σ(⟨φ⟩).

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 volume exclusion interactions only, the free energy calculated from Eq. (1) reduces to an expression determined entirely by the mesostate (the occupancy) N : βH 0 (φ|N ) = -µ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 ⟨φ⟩ takes the following typical sigmoid form as a function of µ:

⟨φ⟩ 0 (µ) = 1/ 1 + e -µ . (7)
The corresponding standard deviation becomes:

σ 0 (µ) = sech (µ/2) 2 √ L = ⟨φ⟩ 0 -⟨φ⟩ 2 0 L . (8)
This latter expression exemplifies the general result that, in the thermodynamic limit, L → ∞, ⟨φ 2 ⟩ → ⟨φ⟩ 2 , and hence, as expected, the standard deviation (of the relative occupancy) tends to zero as L -1/2 .

Fluctuations in the strong cooperativity limit

Away from HF, although ξ 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 below). Taking the limits J → ∞ and µ = 2X -J → -∞ simultaneously, but keeping L finite and X constant, allows Eqs. ( 2) to be simplified to

⟨φ⟩ ∞ (X) = 1 2 1 + tanh L 2ξ ∞ (X) sgn(X) , (9)
and

⟨φ 2 ⟩ ∞ = ⟨φ⟩ ∞ , where ξ ∞ (X) = 1/(2|X|), leading to a large J standard deviation of σ ∞ (X) = 1 2 sech L 2ξ ∞ (X) = ⟨φ⟩ ∞ -⟨φ⟩ 2 ∞ . (10)
The maximal standard deviation at half-filling is due to the large width of the distribution of states with non-zero probability compatible with ⟨φ⟩ = 1/2 (see Section 2.4); in Section 4.2, we show 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.

Fluctuations for intermediate finite cooperativity

The relative mean square occupancy at equilibrium,

⟨φ 2 ⟩ = L -2 L i=1 L j=1 ⟨φ i φ j ⟩ = 1 L 2 Z ∂ 2 µ Z, (11)
is controlled by the 2-point correlation function,

C ij = ⟨φ i φ j ⟩ (with C ii = ⟨φ 2 i ⟩ = ⟨φ i ⟩ = ⟨φ⟩)
. This result can be used as a starting point to gain deeper physical insight into the behavior of the standard deviation σ than the previous result obtained directly from the partition function because it becomes possible to express σ as an explicit function of ⟨φ⟩ and the correlation length ξ.

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 ̸ = j: C ij → ⟨φ i ⟩⟨φ j ⟩ = ⟨φ⟩ 2 . In this J = 0 limit, ⟨φ 2 ⟩ can then be obtained directly from Eq. ( 11) and σ simplifies to σ 0 (Eq. ( 8)). In the other limit of infinite cooperativity (J → ∞) and strong correlations, ξ diverges at HF and away from HF for not too large L remains larger than L/2. Therefore for L finite C ij → ⟨φ 2 i ⟩ = ⟨φ⟩, and σ saturates at σ ∞ (Eq. ( 10)).

Fluctuations in the strong anti-cooperativity limit

By comparing the L = 13 and 14 cases (Figs. 2a and2b), 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 ≪ -1 towards a characteristic limiting shape with global maxima near ⟨φ⟩ = 0.2 and 0.8 and a zero at half-filling (see Fig. 2b). A restricted partition function approach (presented in Section 4.2) that retains only non-overlapping particle-hole pairs leads to a simple but accurate approximation in the limit J → -∞ for L even:

σ even -∞ = L -1/2 ⟨φ⟩ -2⟨φ⟩ 2 for 0 ≤ ⟨φ⟩ ≤ 1/2. ( 12
)
The corresponding standard deviation for 1/2 ≤ ⟨φ⟩ ≤ 1 can be obtained by exploiting particlehole symmetry: σ(⟨φ⟩) = σ(1 -⟨φ⟩). In this limit of strong anti-cooperativity (J → -∞), the correlation length becomes complex, 1/ξ ≃ iπ + 2 cosh(X)e -|J|/2 . Hence, at HF, 1/ξ ≃ iπ, leading directly to an oscillating twopoint correlation function that describes a sequence of non-overlapping particle-hole pairs (antiferromagnetic order in Ising spin language):

C HF i,i+r ≈ 1 + cos(πr) 4 for J → -∞, (13)
where we have taken the real part of C i,i+r and assumed that r/L ≪ 1. C HF i,i+r 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.

Probability distribution function (PDF) of the occupancy and its dependency on ⟨φ⟩ 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 µ, P (N ; J, µ) through Monte Carlo simulations or exact enumeration and study how it depends on the interaction potential J. To study P (N ; J, ⟨φ⟩) for fixed J and ⟨φ⟩, the experimentally accessible quantity, we varied J but kept the mean occupancy fixed (by adapting µ to any given choice of J). Taking L = 13, Fig. 3 shows how the PDF behaves for three different positive values of the interaction potential and three characteristic values of the mean relative occupancy: ⟨φ⟩ = 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, 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 bi-modal form.

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. 3b approximates this Gaussian behavior well because the average value, ⟨N ⟩, is far enough from the boundaries that the tail of the distribution does 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 ⟩ < 2Lσ or L -⟨N ⟩ < 2Lσ (i.e., when the tails of the PDF touch the system boundaries).

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 ∞ (N ; ⟨φ⟩) = (1 -⟨φ⟩)δ N,0 + ⟨φ⟩δ N,L (14)
and the system becomes well described by an effective two-state (bi-modal) system. 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 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. 3a and3c, with an asymmetry that increases with the value of the interaction potential.

PDF in the strong anti-cooperativity limit

The situation is more complicated for strong anti-cooperativity (J ≪ -1) because the system cannot, in general, be reduced to a simple one or two-state system, except near zero and halffilling. Furthermore, near HF, the reduction depends on whether L is even or odd (see Fig. 4). 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. 2b). At HF, there is only one allowed state, consisting of L/2 non-overlapping particle-hole pairs [(φ i , φ 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, σ HF -∞ = 1/(2L), (see PDF for L = 13 in Fig. 4). 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).

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 → ∞), σ = 1/2; (ii) for no cooperativity (J = 0), σ = 1/(2 √ L); (iii) for strong anti-cooperativity (J → -∞) and L odd, σ = 1/(2L); and (iv) for strong anti-cooperativity (J → -∞) and L even, σ = 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 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.

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 [START_REF] Purcell | Life at low reynolds number[END_REF]. They have subsequently evolved a variety of compatible motility mechanisms that have been widely studied [START_REF] Kearns | A field guide to bacterial swarming motility[END_REF][START_REF] Spormann | Gliding motility in bacteria: Insights from studies of myxococcus xanthus[END_REF][START_REF] Henrichsen | Twitching motility[END_REF][START_REF] Martínez | Sliding motility in mycobacteria[END_REF][START_REF] Dvoriashyna | Hydrodynamics and direction change of tumbling bacteria[END_REF][START_REF] Wadhwa | Bacterial motility: machinery and mechanisms[END_REF]. One type of motility involves rotating one or several flagella that propel bacteria through aqueous media [START_REF] Nakamura | Flagella-driven motility of bacteria[END_REF][START_REF] Grognot | More than propellers: how flagella shape bacterial motility behaviors[END_REF][START_REF] Brown | Flagellar hook flexibility is essential for bundle formation in swimming escherichia coli cells[END_REF][START_REF] Gabel | Speed of the flagellar rotary motor of Escherichia coli varies linearly with protonmotive force[END_REF][START_REF] Hosking | Escherichia coli MotAB proton channel unplugged[END_REF][START_REF] Sowa | Hybrid-fuel bacterial flagellar motors in escherichia coli[END_REF][START_REF] Ryu | Torque-generating units of the flagellar motor of escherichia coli have a high duty ratio[END_REF][START_REF] Reid | The maximum number of torque-generating units in the flagellar motor of escherichia coli is at least 11[END_REF][START_REF] Fung | Powering the flagellar motor of escherichia coli with an external voltage source[END_REF]. 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 [START_REF] Berg | The rotary motor of bacterial flagella[END_REF][START_REF] Nirody | The biophysicist's guide to the bacterial flagellar motor[END_REF][START_REF] Fukuoka | Sodium-dependent dynamic assembly of membrane complexes in sodium-driven flagellar motors[END_REF][START_REF] Hosking | Escherichia coli MotAB proton channel unplugged[END_REF][START_REF] Kojima | Mutations targeting the C-terminal domain of FliG can disrupt motor assembly in the Na(+)-driven flagella of Vibrio alginolyticus[END_REF][START_REF] Biquet-Bisquert | The dynamic ion motive force powering the bacterial flagellar motor[END_REF]. 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 [START_REF] Tipping | Load-dependent assembly of the bacterial flagellar motor[END_REF][START_REF] Bai | Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch[END_REF][START_REF] Lele | Mechanism for adaptive remodeling of the bacterial flagellar switch[END_REF][START_REF] Lele | Dynamics of mechanosensing in the bacterial flagellar motor[END_REF]. Its mechanisms of adaptation to external stimuli have been widely studied, becoming a model molecular machine to investigate aspects such as mechano-sensitivity [START_REF] Castillo | The c-terminal periplasmic domain of motb is responsible for load-dependent control of the number of stators of the bacterial flagellar motor[END_REF][START_REF] Bai | Conformational spread as a mechanism for cooperativity in the bacterial flagellar switch[END_REF][START_REF] Lele | Dynamics of mechanosensing in the bacterial flagellar motor[END_REF][START_REF] Nord | Catch bond drives stator mechanosensitivity in the bacterial flagellar motor[END_REF][START_REF] Hu | Molecular stretching modulates mechanosensing pathways[END_REF][START_REF] Wadhwa | Mechanosensitive remodeling of the bacterial flagellar motor is independent of direction of rotation[END_REF], chemotaxis [START_REF] Sourjik | Responding to chemical gradients: bacterial chemotaxis[END_REF][START_REF] Ma | Bacterial flagellar motor switch in response to cheyp regulation and motor structural alterations[END_REF][START_REF] Yuan | Ultrasensitivity of an adaptive bacterial motor[END_REF] and dynamic subunit exchange [START_REF] Tusk | Subunit exchange in protein complexes[END_REF].

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 [START_REF] Santiveri | Structure and function of stator units of the bacterial flagellar motor[END_REF][START_REF] Deme | Structures of the stator complex that drives rotation of the bacterial flagellum[END_REF][START_REF] Zhu | Conformational change in the periplamic region of the flagellar stator coupled with the assembly around the rotor[END_REF][START_REF] Kojima | Stator assembly and activation mechanism of the flagellar motor by the periplasmic region of motb[END_REF][START_REF] Tipping | Quantification of flagellar motor stator dynamics through in vivo proton-motive force control[END_REF][START_REF] Kojima | The helix rearrangement in the periplasmic domain of the flagellar stator b subunit activates peptidoglycan binding and ion influx[END_REF][START_REF] Terahara | Na-induced structural transition of motps for stator assembly of the flagellar motor[END_REF][START_REF] Lin | Flil association with flagellar stator in the sodium-driven vibrio motor characterized by the fluorescent microscopy[END_REF]. 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 ≈ 13 stator units and that the system is mechanosensitive in that stator unit stoichiometry scales with the external torque induced by the viscous drag upon the flagellum turning inside an aqueous medium [START_REF] Reid | The maximum number of torque-generating units in the flagellar motor of escherichia coli is at least 11[END_REF].

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, ω, and the torque produced, τ (from the relation τ = γω, where γ 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, ⟨N ⟩, by varying the beads' size and hence the viscous load, because the binding of stators to the BFM is mechanosensitive (see, e.g., [START_REF] Nord | Catch bond drives stator mechanosensitivity in the bacterial flagellar motor[END_REF], 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 [START_REF] Nord | Catch bond drives stator mechanosensitivity in the bacterial flagellar motor[END_REF][START_REF] Perez-Carrasco | Relaxation time asymmetry in stator dynamics of the bacterial flagellar motor[END_REF]).

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. 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, µ 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, ε, depends on the viscous load (in our case the size of the bead and the viscosity of the surrounding medium). Load-induced changes in ε 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. 5 shows a comparison between the theoretical standard deviation curves, already presented in Fig. 2a, 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 certain level of cooperativity in the system, J ≈ 1.5, is consistent with the experimental observations, 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. 6 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 we don't observe a clear 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 the experiment. 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 bi-modal (with only the completely full and empty states having significant probabilities) for slightly higher values of the interaction potential, J ≳ 5, behavior that is not seen in the experimental PDF data. We see from Section 2.4 (Fig. 3) that for the chosen system size (L = 13), J ≈ 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.

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, 2 k B T (i.e., J ≈ 2), for the short-range 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 model, using the mapping (59) derived in Section 4.2, the corresponding value is

J IR = 2J SR /(L -1) ≈ 0.33].
For the characteristic size of the BFM (≈ 10 binding sites), slightly higher values of cooperativity would lead to a motor that is bi-modal in the occupancy and, therefore, would exhibit switch-like behavior for the produced torque not compatible with the required motor characteristics.

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 ε 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 µ and therefore a modification of the average occupancy and standard deviation (fluctuations).

For the short-range model, these modifications to µ (extracted from the experimental data) are shown in Fig. 7 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 µ window in this way we observe that, for 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 [START_REF] Perez-Carrasco | Relaxation time asymmetry in stator dynamics of the bacterial flagellar motor[END_REF]); 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).

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 (J ≈ 2) are optimal given the µ 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. 7, we observe that if a two-state switch-like motor operation were sought after in order to cover a wider range of occupancy, then higher values of J would be the best choice, the price to pay would be a strong increase in fluctuations. 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 was 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. One important question concerns how to integrate cooperativity into kinetic models that allow one to account for the relaxation time asymmetry between stall and resurrection [START_REF] Perez-Carrasco | Relaxation time asymmetry in stator dynamics of the bacterial flagellar motor[END_REF].

Figure 7: Values of the effective chemical potential µ extracted from the experimental data for different values of the interaction potential J (L = 13). We situate them in the contour plots of ⟨φ⟩ and σ in order to find the interval of chemical potential in which the BFM would work for fixed J.

Methods

1D short-range lattice gas (SRLG)

Analytical solutions for the equilibrium state can be obtained using the transfer matrices formalism. According to it, the partition function of the system can be written as

Z = {φ i } exp[-βH(φ)] = {φ i } L i=1 ⟨φ i | T |φ i+1 ⟩ , (15)
where T is the transfer matrix, which for the system we are working on, is

T = 1 e µ 2 e µ 2 e (J+µ) . (16)
In the case of periodic boundary conditions, equation ( 15) can be simplified to Z = Tr T L = λ L + + λ L -, being λ + and λ -the eigenvalues of the transfer matrix, given by:

λ ± = e X (cosh X ± sinh 2 X + e -J ), (17)
where X = β 2 (J +µ). From this point, we can obtain the mean occupancy and mean square values by continuing to use the transfer matrix formalism or following a thermodynamic approach.

According to the grand-canonical formalism, the average occupancy number of a system is given by the expression

⟨φ⟩ = - 1 L 2 ∂ Ω ∂µ , ( 18
)
where Ω is the grand potential, which depends on the partition function as Ω = -L β ln Z. Let's remind that Z = λ L + + λ L -and that λ ± dependency on J and µ is given by equation [START_REF] Godreche | Response of non-equilibrium systems at criticality: ferromagnetic models in dimension two and above[END_REF]. Using the chain rule, we can find that

⟨φ⟩ = 1 βZ λ (L-1) + ∂ λ + ∂X ∂ X ∂µ + λ (L-1) - ∂ λ - ∂X ∂ X ∂µ , (19) being
∂ X ∂µ = β 2 and ∂ λ ± ∂X = λ ± (1 ± sinh X √ sinh 2 X+e -βJ ).
Replacing these in the previous equation, we find the expression

⟨φ⟩ = 1 2 1 + tanh L 2ξ sinh X sinh 2 X + e -βJ , (20)
where ξ is the correlation length defined previously. As clarification, tanh

( L 2ξ ) = λ L + -λ L - λ L + +λ L - .
Once the calculation of the mean occupancy, we may calculate the mean square value, given by

⟨φ 2 ⟩ = 1 L 2 1 β 2 Z ∂ 2 Z ∂µ 2 . (21)
By carrying out this calculation, we obtain the expression

⟨φ 2 ⟩ = 1 4 1 + sinh 2 X sinh 2 X + e -J + tanh L 2ξ 2 sinh X sinh 2 X + e -J + 1 L e -J cosh X (sinh 2 X + e -J ) 3/2 . (22)
This way, we can obtain the analytical expression for the standard deviation, which is σ = ⟨φ 2 ⟩ -⟨φ⟩ 2 . The regime of thermodynamic limit corresponds to L → ∞, for which tanh (L/2ξ) → 1, so we may find the expressions on the thermodynamic limit as

⟨φ⟩ = 1 2 1 + sinh X sinh 2 X + e -J , (23)
and

⟨φ 2 ⟩ = 1 4 1 + sinh X sinh 2 X + e -J 2 . ( 24
)
That way, as in the case with the equations of the transfer matrix method, we have that ⟨φ 2 ⟩ = ⟨φ⟩ 2 . Hence the standard deviation on the thermodynamic limit, σ, is zero.

Probability distribution function of the occupancy

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

P ∞ (N ; ⟨φ⟩) = (1 -⟨φ⟩)δ N,0 + ⟨φ⟩δ N,L (25)
and the system becomes well described by an effective two-state system. The standard deviation for the occupancy N therefore saturates at Lσ ∞ (= L/2 at half-filling), where we recall that σ ∞ = ⟨φ⟩ -⟨φ⟩ 2 , and it becomes easy to calculate all moments of N : ⟨N m ⟩ = L m ⟨φ⟩ (or ⟨φ m ⟩ = ⟨φ⟩). The moments are defined by

µ n = ⟨(φ -⟨φ⟩) n ⟩ (26)
and the standard moments by µ n /σ n , where σ = √ µ 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 γ = µ 3 /σ 3 ,

γ ∞ = µ 3,∞ σ 3 ∞ = ⟨φ⟩ -3⟨φ⟩ 2 + 2⟨φ⟩ 3 σ 3 ∞ (27)
and kurtosis, κ = µ 4 /σ 4 ,

κ ∞ = µ 4,∞ σ 4 ∞ = ⟨φ⟩ -4⟨φ⟩ 2 + 6⟨φ⟩ 3 -3⟨φ⟩ 4 σ 4 ∞ . (28)
These results for the effective two-state system are very different from those predicted for a Gaussian PDF.

For strong anti-cooperativity (J ≪ -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. 4). 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 -∞ (N ; 1/2) = δ N,L/2 (L even), (29)
(see Fig. 4) 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 -∞ (N ; 1/2) = 1 2 δ N,(L-1)/2 + 1 2 δ N,(L+1)/2 (L odd), (30)
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 (nonparticle-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 partition function approach consists in keeping in the sum over all possible microstates {φ 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[-βH 0 (N )] = e µN multiplied by a multiplicity (or binomial coefficient) C L N = L!/ [N ! (L -N )!] (related to the configurational entropy) that gives the number of microstates φ consistent with occupancy N :

Z 0 = {φ i } e -βH 0 {φ i } = L N =0 C L N e µN = (1 + e µ ) L . (31)
From Z 0 we immediately obtain the expected J = 0 results:

⟨φ⟩ 0 = (LZ 0 ) -1 ∂ µ Z 0 = 1/ 1 + e -µ (32)
and

σ 0 = L -1 ∂ µ ⟨φ⟩ 0 = L -1/2 ⟨φ⟩ -⟨φ⟩ 2 . ( 33
)
By introducing the usual (configurational) entropy of mixing,

S mix = -k B L[φ ln(φ) + (1 -φ) ln(1 -φ)], (34)
the multiplicity C L N can be accurately approximated (using Stirling's formula) for L ≫ 1 and 0 < φ = N/L < 1 by

C L N ≈ e Ls mix (φ) 2πLφ(1 -φ) , (35)
where

s mix ≡ S mix /(Lk B ) = -[φ ln(φ) + (1 -φ) ln(1 -φ)] (36)
is a dimensionless entropy of mixing. The entropy of mixing is a non-monotonic (concave) function of φ that goes to 0 for zero and full filling and reaches a maximum at HF. In Eq. ( 35) for L ≫ 1, the main variation comes from the exponential, and it is possible under certain conditions to replace φ in the prefactor by ⟨φ⟩ without loss of accuracy. The PDF for J = 0 can then be accurately approximated by

P 0 (N ; ⟨φ⟩) ≡ C L N e µN Z 0 ≈ e -βF 0 (φ;⟨φ⟩) (2π) 1/2 σ 0 LZ 0 (⟨φ⟩) , (37)
where

βF 0 (φ; ⟨φ⟩) = -L φ ln ⟨φ⟩ 1 -⟨φ⟩ + s mix (φ) (38)
is an effective free energy for φ and we have written

Z 0 (⟨φ⟩) = 1 (1 -⟨φ⟩) L , (39)
and

µ(⟨φ⟩) = ln ⟨φ⟩ 1 -⟨φ⟩ (40)
as functions of ⟨φ⟩.

For L ≫ 1, we can use the saddle point approximation to obtain a simple Gaussian form for P 0 (N ; ⟨φ⟩): by expanding βF 0 (φ; ⟨φ⟩) around φ = ⟨φ⟩ (the position of its minimum) to second order in (φ -⟨φ⟩) we find

P 0 (N ; ⟨φ⟩) ≈ exp -(φ-⟨φ⟩) 2 2σ 2 0 (2π) 1/2 σ 0 L , (41)
where the variance

σ 2 0 = -[Ls ′′ mix (⟨φ⟩)] -1 = (⟨φ⟩ -⟨φ⟩ 2 )/L ( 42
)
is related to the second derivative of the entropy of mixing, or curvature, evaluated at ⟨φ⟩. Since -s ′′ mix (⟨φ⟩) has a minimum at half-filling, σ 0 has a maximum there.

The Gaussian approximation for P 0 (N ; ⟨φ⟩) 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 φ from -∞ to +∞ about its maximum value at ⟨φ⟩. The Gaussian approximation is itself clearly only valid when the tails of the distribution are far enough from the extremal values (φ = 0 and 1), which is the case for the three J = 0 PDFs displayed in Fig. 3 (for ⟨φ⟩ = 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 L N , Eq. ( 35), should be accurate over the whole range of φ except close to the extremal values.

In the strong cooperativity limit (J → +∞), only the empty and full states need be retained in the restricted partition function, Z +∞ , 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

Z +∞ (X; µ) = 1 + e 2XL , (43)
where X = (J + µ)/2 is considered to be kept fixed as J → +∞ and µ → -∞. From Z +∞ (X; µ) we obtain

⟨N ⟩ +∞ = Z -1 +∞ ∂ µ Z +∞ = Le 2XL /Z +∞ , (44)
which allows us to read off the expected results,

P +∞ (0; ⟨φ⟩) = Z -1 +∞ = 1 -⟨φ⟩ (45)
and

P +∞ (L; ⟨φ⟩) = e 2XL Z -1 +∞ = ⟨φ⟩, (46)
directly from Z +∞ [within this approximation the PDF P +∞ (N ; ⟨φ⟩) vanishes for all other values of N ].

For strong anti-cooperativity (J ≪ -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 ≫ 1) limit where only two states need be retained.

We start with the simpler L even case for which the system can be approximated by nonoverlapping, 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 partition function for 0 ≤ ⟨φ⟩ ≤ 1/2 in this particle-hole pair approximation is

Z even -∞ = L/2 N d =0 D L N d e µN d , (47)
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 L N d counts the number of distinct ways of putting N d pairs on a periodic lattice of size L. Since calculating D L N d 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 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 partition function

Z even d = L/2 N =0 2C L/2 N d e µN d = 2 (1 + e µ ) L/2 , (48)
where

C L/2
N d 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. ( 47) to Eq. ( 48), we have replaced

D L N d by 2C L/2
N d and therefore simplified the problem to that of two independent non-interacting quasi-particle (i.e., particle-hole pair) systems [cf. Eq. ( 31))]. Following the above discussion for the J = 0 case, we can immediately conclude that ⟨φ⟩ = ⟨N d /L⟩ ≈ ⟨φ⟩ even d , where

⟨φ⟩ even d = (LZ d ) -1 ∂ µ Z even d = 1/[2 1 + e -µ ] (49)
and therefore σ even -∞ ≈ σ even d , where

σ even d = L -1 ∂ µ ⟨φ⟩ even d = L -1/2 ⟨φ⟩ -2⟨φ⟩ 2 . ( 50
)
The replacement of D L N d by 2C

L/2 N d 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 Z even d 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 Z even d 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.

Inspection of Fig. 2b shows that the approximate form σ even d 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 σ even d , however, the approximate result cannot capture the lack of symmetry about 1/4 filling: the observed skewness towards lower fillings seen in Fig. 2b, and more clearly seen in the PDFs for L = 14 and ⟨φ⟩ = 1/3 presented in Fig. 4, arises because for 1 < N d < L/2 the approximation 2C L/2 N d underestimates the exact particle-hole pair multiplicity D L N d more severely for lower fillings than for higher ones, because 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. 4 for ⟨φ⟩ = 1/3 can be approximated using the approximate particle-hole pair approach that maps the problem onto free quasi-particles (pairs). In this case we can adapt the results of the J = 0 case to estimate the particle-hole pair PDF for J → -∞ and L even,

P even d (N ; ⟨φ⟩) ≡ C L/2 N d e µ even d N d Z even d , (51)
using the Gaussian approximation:

P even d (N ; ⟨φ⟩) ≈ exp -(φ-⟨φ⟩) 2 2(σ even d ) 2 (2π) 1/2 σ even d L , (52)
where we have used

Z even d (⟨φ⟩) = 2 (1 -2⟨φ⟩) L/2 , (53)
and

µ even d (⟨φ⟩) = ln 2⟨φ⟩ 1 -2⟨φ⟩ . (54)
The Gaussian PDF approximation is very accurate compared with P even d (N ; ⟨φ⟩) for L = 14 and ⟨φ⟩ = 1/3, although P even d (N ; ⟨φ⟩) itself underestimates the exact maximum by about 30% and cannot capture the skewness observed in the exact results (see Fig. 4). The particle-hole pair approximation does reproduce accurately, however, the width of the exact PDF.

In the strong anti-cooperativity limit (J → -∞) 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 partition function for 0 ≤ ⟨φ⟩ ≤ 1/2 in this particle-hole pair approximation is

Z odd -∞ = (L-3)/2 N d =0 O L N d e µN d + Le µ(L-1)/2 + Le -|J|+µ(L+1)/2 , (55)
where (L -1)/2 is the maximum number of non-overlapping pairs on a lattice of odd size L and the multiplicity O L N d 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 L N d 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 ≤ N d ≤ (L -3)/2, i.e., those that have pairs on both lattices, we obtain the following approximate restricted particle-hole pair partition function

Z odd d = (L-3)/2 N d =0 2C (L+1)/2 N d e µN d + Le µ(L-1)/2 + Le -|J|+µ(L+1)/2 , (56)
which can be evaluated by completing the binomial sum and then subtracting the added terms:

Z odd d = 2 (1 + e µ ) (L+1)/2 -e µ(L-1)/2 -e µ(L+1)/2 (2 -Le -|J| ) (57)
where C , an approximate results that could then be compared with the exact ones. We will not, however, pursue this approach any further here.

Infinite range (IR) model

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

Z = {φ i } e -H 0 -H int = Z 0 ⟨e -H int ⟩ 0 ≈ Z 0 e -⟨H int ⟩ 0 , (58)
where the last line gives the first order cumulant expansion. The subscript 0 indicates a statistical average with respect to the non-interacting (Langmuir) model (ideal Lattice Gas). By calculating ⟨H int ⟩ 0 , we obtain ⟨H int ⟩ 0 = -J SR L⟨φ⟩ 2 0 for the SR model and -1 2 J IR L(L -1)⟨φ⟩ 2 0 for the IR model, results that suggest the following mapping between the two models:

J SR ↔ 1 2 J IR (L -1). ( 59
)
This mapping is exact to lowest order in the coupling constant. 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 (Langmuir) model. This result can be extended to find we find

P I R(N ; ⟨φ⟩) ≈ exp -(φ-⟨φ⟩) 2 2σ 2 IR (2π) 1/2 σ IR L , (60)
where the variance σ IR can be approximated using the approximation for the SR model and the mapping.

We adopt the following strategy to avoid carrying out an explicit calculation of the standard deviation σ 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 ⟨φ⟩(µ) and the variance v(µ) = σ 2 (µ) in powers of J. To first order in J SR we find

σ 0 σ SR 2 = 1 -2J SR Lσ 2 0 + O(J 2 SR ) = 1 -2J SR σ 2 ∞ + O(J 2 SR ). ( 61
)
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:

σ 0 σ IR 2 = 1 -L(L -1)J IR σ 2 0 + O(J 2 IR ) = 1 -(L -1)J IR σ 2 ∞ + O(J 2 IR ). ( 62
)
One can check the above method for SR model at half-filling:

1 - σ 0 σ SR 2 = 1 2 J SR + O(J 2 SR ). (half-filling) (63)
It would be useful to find a good interpolation formula for σ as an explicit function of L, J, and ⟨φ⟩ 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,

σ 0 σ SR 2 = exp(-2J SR σ 2 ∞ ) × 1 + O( σ 2 ∞ J SR L ) (half-filling) (64)
(where σ 2 ∞ = 1/4 at HF), which shows that it's best to use an exponential resummation of the small J perturbation expansion.

Experimental data

Figure 1 :

 1 Figure 1: Contour plots of the relative mean occupancy in equilibrium ⟨φ⟩ (left), and the standard deviation of the relative occupancy in equilibrium σ = ⟨φ 2 ⟩ -⟨φ⟩ 2 (right), as functions of the dimensionless interaction potential J and chemical potential µ, according to expressions from Eq. (2) and (4), for a system of size L = 13.

Figure 2 :

 2 Figure 2: Parametric plots of the standard deviation of the relative occupancy, σ, versus the mean relative occupancy, ⟨φ⟩, 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.

Figure 3 :

 3 Figure 3: 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 an exact counting 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.

Figure 4 :

 4 Figure 4: Discrete equilibrium probability distribution of the number of bound particles, P (N ), for J = -10 and three different mean occupancies: ⟨φ⟩ = 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).

Figure 5 :

 5 Figure 5: Parametric plot of the standard deviation of the relative occupancy, σ, versus the mean occupancy, ⟨φ⟩, comparing the experimental values with the theory; here, L = 13. The experimental data corresponds to five different applied viscous loads with three different microparticles (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).

Figure 6 :

 6 Figure 6: 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.

  L+1)/2 N d counts the number of distinct ways of placing N d particles on a lattice of size (L + 1)/2. Z odd d correctly counts the multiplicity of the (L ± 1)/2 states by construction. Although it 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, Z odd d , can be used to calculate ⟨φ⟩ odd d = (LZ odd d ) -1 ∂ µ Z odd d and σ odd d = L -1 ∂ µ ⟨φ⟩ odd d

Figure 8 :

 8 Figure 8: Contour plot of the relative mean occupancy in equilibrium, ⟨φ⟩, as a function of the dimensionless interaction potential J and chemical potential µ for a system of size L = 14.

Figure 9 :

 9 Figure 9: Contour plot of the standard deviation of the relative occupancy at equilibrium, σ = ⟨φ 2 ⟩ -⟨φ⟩ 2 , as a function of the dimensionless interaction potential J and chemical potential µ for a system of size L = 14.

Figure 10 :

 10 Figure 10: Standard deviation of the mean occupancy at half-filling, ⟨φ⟩ = 1/2, as a function of the dimensionless interaction potential J for a system of size L = 13 (blue) and L = 14 (orange).

Figure 11 :

 11 Figure 11: Comparison between the exact result for the standard deviation (solid lines) and the approximate expression that depends on ⟨φ⟩, Eq. (66) (dashed lines).

Table 1 :

 1 Data of average occupation and standard deviation corresponding to Fig.5. Data are obtained by averaging over several traces. As there is not enough traces for good statistics (from 30 to 40), data is also averaged over several time points at equilibrium.

	Bead size (drag)	⟨φ⟩	σ
	300 nm	0.29 ± 0.03	0.17 ± 0.02
		0.355 ± 0.006 0.153 ± 0.011
		0.41 ± 0.02	0.15 ± 0.02
	300 nm + Glycerol 0.461 ± 0.005 0.175 ± 0.007
	500 nm	0.56 ± 0.02	0.18 ± 0.02
		0.592 ± 0.011	0.21 ± 0.02
		0.67 ± 0.02	0.18 ± 0.04
	500 nm + Glycerol 0.705 ± 0.005 0.158 ± 0.005
	1300 nm	0.74 ± 0.03	0.13 ± 0.02
		0.80 ± 0.03	0.153 ± 0.012
		0.828 ± 0.006	0.18 ± 0.04
		0.826 ± 0.011 0.161 ± 0.009

Table 2 :

 2 Data of probability distribution corresponding to Fig.6. Data are obtained by averaging over several time points at equilibrium the state of all traces.

			Bead size
		300 nm	500 nm	1300 nm
	# Stator units		Probability
	0	0.12 ± 0.04	0.005 ± 0.014	0.01 ± 0.20
	1	0.09 ± 0.04	0.005 ± 0.013	0
	2	0.11 ± 0.05	0	0
	3	0.15 ± 0.04	0.01 ± 0.02	0
	4	0.15 ± 0.04	0.05 ± 0.02	0
	5	0.22 ± 0.04	0.04 ± 0.03	0
	6	0.13 ± 0.03	0.17 ± 0.05	0
	7	0.033 ± 0.014	0.29 ± 0.10	0.002 ± 0.011
	8	0	0.13 ± 0.04	0.12 ± 0.04
	9	0	0.13 ± 0.04	0.14 ± 0.03
	10	0	0.05 ± 0.02	0.28 ± 0.07
	11	0	0.04 ± 0.02	0.25 ± 0.06
	12	0	0.06 ± 0.02	0.18 ± 0.05
	13	0	0.03 ± 0.02	0.02 ± 0.03

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.

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),

to get a smooth interpolation. The above 5/2 root mean square average leads to our best approximation for σ SR :

This is a convenient form because the dependence on ⟨φ⟩ (through σ ∞ ), J SR , and L is made explicit (σ 0 = L -1/2 σ ∞ ). The crossover to the large J saturation limit is centered on

J co depends on L and filling ⟨φ⟩, increasing slowly with L and reaching a minimum at half-filling for fixed L. The convergence of σ as a function of ⟨φ⟩ to the large J saturation limit is uniform only for J < J co . This approximation is very accurate for all fillings (except possibly very close to zero and full-filling) as long as L -5/2 ≪ 1 (the case for L ≈ 10) and J < J co (see Fig. 11). For L -5/2 ≪ 1 and J > J co , this approximation is accurate only within a window centered on half-filling. Because of the factor σ 2 ∞ in the argument of the exponential, convergence to σ ∞ is fastest at HF and becomes slower and slower as moves away from HF to zero and full filling. For L ≈ 10, the window of accuracy for all J is 0.2 < ⟨φ⟩ < 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 σ growing exponentially with J until the system reaches the crossover region, J ∼ J co , before saturating at σ ∞ for J ≫ 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 (59). 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 (Langmuir) model. This result can be used to find an extended Gaussian approximation for the IR model,

where the variance σ IR can be approximated using the approximation for the SR model (66) and the mapping (59).

Supplemental plots