Sensory adaptation in a continuum model of bacterial chemotaxis—working range, cost-accuracy relation, and coupled systems

Sensory adaptation enables organisms to adjust their perception in a changing environment. A paradigm is bacterial chemotaxis, where the output activity of chemoreceptors is adapted to different baseline concentrations via receptor methylation. The range of internal receptor states limits the stimulus magnitude to which these systems can adapt. Here, we employ a highly idealized, Langevin-equation based model to study how the finite range of state variables affects the adaptation accuracy and the energy dissipation in individual and coupled systems. Maintaining an adaptive state requires constant energy dissipation. We show that the steady-state dissipation rate increases approximately linearly with the adaptation accuracy for varying stimulus magnitudes in the so-called perfect adaptation limit. This result complements the well-known logarithmic cost-accuracy relationship for varying chemical driving. Next, we study linearly coupled pairs of sensory units. We find that the interaction reduces the dissipation rate per unit and affects the overall cost-accuracy relationship. A coupling of the slow methylation variables results in a better accuracy than a coupling of activities. Overall, the findings highlight the significance of both the working range and collective operation mode as crucial design factors that impact the accuracy and energy expenditure of molecular adaptation networks.


I. INTRODUCTION
The surroundings of all living organisms constantly change and diverse mechanical, electrical, and chemical signaling pathways enable organisms to continuously detect these changes [1][2][3][4].Reliable detection of changes in different signal backgrounds requires the adjustment of the sensitivity, which is called sensory adaptation [5].A paradigm of sensory adaptation is chemotaxis of the bacterium Escherichia coli (E. coli ) [6].Chemotaxis allows E. coli to navigate up or down chemical concentration gradients, e.g., to find an environment with higher sugar concentrations.Chemotactic movement is a result of a stochastic control of the rotation direction of the flagella that drive bacterial swimming [3].E. coli swims in a pattern that comprises of successive straight runs and tumble phases, during which it reorients.As the bacterium moves towards higher concentrations of a chemoattractant, the tumbling frequency is reduced, which leads to increasingly long runs that take the cell in the direction of higher concentration.Conversely, the tumbling frequency is increased as the bacterium swims away from a repellent.The presence of ligands is sensed by transmembrane receptors that are called methyl-accepting chemotaxis proteins.Depending on ligand binding, receptors can either be in an active or an inactive state.Chemoattractant binding inhibits receptor activity, which suppresses tumbling by lowering the phosphorylation state of a response regulator CheY [7].Simultaneously, an adaptation of the receptors to the ligand concentration occurs via methylation or demethylation of their four to five glutamate residues with the help of the methyltransferase CheR and the methylesterase CheB, respectively [8,9].CheR pref-erentially methylates inactive receptors and increases receptor activity through methylation if receptor activity is reduced by an attractant.CheB preferentially demethylates active receptors and thereby lowers activity upon removal of attractants.Consequently, adaptation restores the receptor activities to their original value after an initial response to a ligand concentration change [5].This feedback mechanism can in principle ensure a precise and robust detection of concentration gradients, regardless of the ambient background concentration of ligands.However, experiments also revealed a limited working range of adaptation with imprecise adaptation to very high concentrations of attractants [10][11][12].This limitation most likely results from a saturation of available methylation sites [13].While the impact of a limited working range on the sensory performance have been studied theoretically [11,13], it is less well-understood how the energetic cost of adaptation depends on the working range of the mechanisms.Chemosensory adaptation that is realized by a finite number of receptor methylation sites can be idealized as a dynamical process with a bounded range of internal variables.How such a bounded internal phase space affects the energetics of nonlinear biological processes is generally an open question.This issue motivates our study of sensory adaptation with a highly idealized model where the interplay of phase-space boundaries of internal variables, energetics, and adaptive performance can be studied.Chemotaxis receptor molecules form clusters at the cell membrane [14,15], where individual molecules form homodimers and three such homodimers form a trimer [16][17][18].A receptor core complex, the minimal unit of signaling, contains two trimers and shows a non-linear doseresponse [19].Assemblies of multiple chemoreceptors work in a cooperative manner to enhance the sensitivity of the response to subtle concentration gradients [20][21][22].On a molecular level, cooperativity can result from enzymes that change the methylation state of the dimer they are bound to as well as the methylation state of neighboring dimers, e.g., CheR and CheB [23,24].The molecular mechanisms of receptor cooperativity, as well as its functional advantages, e.g., to achieve optimal signal-to-noise ratios, have been studied extensively.However, the impact of receptor cooperativity on the energy consumption required for adaptation is hardly explored and warrants further theoretical research.
A variety of models for sensory adaptation in chemotaxis have been proposed over the past decades [25][26][27][28][29][30][31].Idealized feedback systems can be represented by only a few coupled differential equations [32,33].However, chemosensory adaptation is usually modeled as a stochastic process since the number of molecules constituting a sensory apparatus in individual cells is typically quite small and thermal noise and active processes result in noisy signals.This noise affects the reliability of information transmission between an input signal and the output [34][35][36] and thereby determines the precision of the cellular response.However, fluctuations are also intimately connected with the energetic cost of adaptation.From an information-theoretic viewpoint, sensory adaptation always entails a thermodynamic cost since it implies an irreversible erasure of stored information along with a measurement [37].Moreover, the two fundamental models of noisy cellular adaptation, namely the negative feedback mechanism in chemotaxis and the incoherent type-1 feed-forward mechanism, only show adaptation if the systems' states are held in non-equilibrium [38,39].The energetic cost of maintaining a non-equilibrium state required for adaptation can be related to sensory performance by an energy-speed-accuracy relation [38], which states that the dissipation rate is proportional to the negative logarithm of the adaptation error, both measured in steady state after application of a constant stimulus.
Most work on the energy-speed-accuracy relationship is based on a discrete description of receptor states using master equations [38], for which the original relationship can be generalized [40].Here, one assumes two discrete states for the fast receptor activity and four to five discrete states that represent the slowly changing methylation level that acts as a a negative control element in E. coli chemotaxis.Alternatively, the receptor states can be represented in a continuum picture using Langevin equations for both the response and the control variables.Based on simulation results for such a model, the generality of the energy-speed-accuracy relation has been questioned [41].The Langevin-equation based model is highly idealized and not as directly applicable to chemotaxis as the discrete model.However, it has the advantage that the adaptation process and the associated fluxes of probability can be clearly visualized in a two-dimensional phase space.This facilitates an intuitive understanding of the adaptation process.In principle, one may even study a hybrid model where the states of the control variable of the receptor are discrete and the response variable is continuous, although the technical advantage of such an approach is not obvious.
In this article, we take advantage of the clear representation offered by the Langevin-equation approach to study the relationships between adaptation, working range, and energetic cost for individual and coupled systems.Unlike earlier studies [38,40], we directly simulate the Langevin dynamics, taking into account the phasespace boundaries.Plotting the adaptation error against the dissipation rate for varying input signals, we find that the finite working range results in a linear cost-accuracy relation, which highlights the importance of the finite working range of the internal variables on the performance and the cost of adaptation.This result complements the well-known logarithmic energy-speed-accuracy relation [38].Also, to underpin our numerical findings, we derive an analytical expression for the energetic cost of adaptation as a function of input stimulus strength.Furthermore, motivated by the discovery that receptor clusters consist of minimal sensory units with a non-linear response [19], we consider pairs of receptors in the spirit of Ref. [42] and ask how a linear coupling of the sensor states affects the overall cost-accuracy relationship.To this end, we study two examples of coupled systems and find that, in general, a linear coupling reduces the cost of adaptation without affecting its accuracy.

A. Model Introduction
An idealized model of bacterial adaptation can be built from a three-node network with a negative feedback.We focus on a Langevin-type model suggested in Ref. [38], where all variables are assumed to have non-negative, real values.The model comprises of a variable a for the rapidly changing receptor activity and a slowly changing control variable m, representing the receptor methylation state.The activity a responds to changes in an input stimulus s and also represents the output of the system.
Following a stimulus, a change in the activity a drives a slow compensatory change in the control variable m, which ultimately brings a back to the stimulusindependent setpoint a 0 , see FIGs.1a and 1b.The dynamics of the activity and the control (methylation) variables is described by a pair of coupled Langevin equations as Here, η a (t) and η m (t) both represent white noise.Denoting expectation values as ⟨. ..⟩, the noise satisfies ⟨η i (t)η j (t ′ )⟩ = 2∆ i δ ij δ(t ′ − t) where the indices i and j represent a or m.The functions F a and F m provide a lumped description of the nonlinear biochemical reac-tions that govern the dynamics of fast response and slow adaptation, respectively [38].
The Fokker-Planck equation (FPE) that corresponds to the Langevin equations (1) describes the evolution of the phase-space probability distribution function (PDF) P ≡ P (a, m, t) as where J a ≡ F a P − ∆ a ∂ a P and J m ≡ F m P − ∆ m ∂ m P are the components of the probability current J.
It can be shown that the system must obey the following condition to satisfy detailed-balance which is also derived in the Supplementary Material.This detailed-balance condition for the negative feedback control mechanism implies that adaptation cannot be realized in equilibrium and is inherently dissipative [43].
We specify the functional form of the forces driving the dynamics of a and m as and where is assumed to be a constant and β ∈ [0, 1] is a parameter that determines whether the system is in equilibrium, β = 0, or out-of-equilibrium, β > 0. The rate of the response and adaptation processes are determined in Eqs.(4,5) by the constants ω a and ω m , respectively.We assume that such that the time scales of response and adaptation are well-separated.
The function G(s, m) is chosen to be a Hill equation in the Michaelis-Menten limit, i.e., with Hill coefficient where we set K 0 = 1.The functional form chosen in Eq. (8) ensures that ∂ m G > 0 and ∂ s G < 0, which is required for the negative feedback mechanism [38].In a deterministic setting, for any β, there exists a fixed point where the two nullclines of the deterministic system intersect.Consequently, there is a value of the control variable m = m * that satisfies Moreover, there exists a β c such that the fixed point is stable only if β > β c , where 0 < β c < 1 [38].Thus, stable adaptation is only achieved if β > β c .Phase portraits of the system with stable fixed points for β = 1 and unstable fixed points for β = 0 are shown in FIGs.1c and 1e, respectively.The corresponding probability distributions are presented in FIGs.1d and 1f.Throughout this article, we focus for simplicity on the case of β = 1, which is called the fully adaptative state which represents the perfect adaptation limit.

Phase space & Adaptation Error
For the model described in Eqs.(1), the activity a and control state m are stochastic variables that lie in a phase space denoted as Ω ≡ [0, 1]×[0, m 0 ], with m 0 corresponding, e.g., to the maximum number of methylation sites in a receptor.To solve the FPE, Eq. ( 2), on the interior of the phase space, Ω ≡ Ω \ ∂Ω, a set of well-defined boundary conditions is required.For this model, we impose reflective boundary conditions, which state that the current vector obeys where n is the outward-pointing vector normal to the boundary ∂Ω.
The setpoint for the activity a 0 is assumed to be a constant.Since the function F m is designed to be proportional to (a − a 0 ), the generalized force driving methylation changes vanishes at the setpoint.The integral feedback mechanism ensures that the statistical average of the activity ⟨a(t)⟩ = Ω a P (a, m, t) da dm (11) approaches a 0 in the steady state [38].Denoting the steady-state average as ⟨. ..⟩ ss , the adaptation error quantifies the accuracy of adaptation in the stationary limit after the decay of transient responses to a signal.

Effective Temperature & Energy Dissipation Rate
In equilibrium, the fluctuation-dissipation theorem (FDT) can be used to relate the temperature of the system to the noise amplitude and the damping in the system [44].In the following, we set the Boltzmann constant to unity, k B ≡ 1, and define the effective temperatures pertaining to the stochastic variables of the system in Eq. ( 1) as follows We assume that the response and adaptation processes described in Eq. ( 1) are coupled to the same heat bath.Therefore, the effective temperatures are equal T eff(m) = T eff(a) = T , which entails that C = 1.
Since the product of temperature and the entropy production rate equals the dissipation rate in a stationary, isothermal setting, see Ref. [45], the stationary dissipation rate can be expressed as This equation for the dissipation rate is only valid for a system with single heat bath and temperature T .For a system with multiple heat baths, additional terms must be included in Eq. ( 14) to account for the heat flux between reservoirs [46,47].An analytical approximation for the dissipation rate was given in Ref. [38], where Eq. ( 7) and Laplace's integral method were used to estimate the integrals in Eq. ( 14), leading to the following expression where C = 1, the superscript "a" stands for the analytical approximation proposed in Ref. [38], and m * is the mcoordinate of the fixed point of the deterministic system i.e., it satisfies Eq. ( 9).Also, employing the framework of stochastic thermodynamics [45,46], the dissipation rate for an individual trajectory "t" of the system can be expressed as where • indicates a product in the Stratonovich sense [46].For an ensemble of trajectories, we take the statistical average of Eq. ( 16) to obtain an expression for the dissipation rate in stationary state as In contrast to Eq. ( 14), Eq. ( 17) can be straight-forwardly evaluated using the trajectories generated by simulations of the stochastic dynamics.

C. Simulation Methods
To simulate the system of Langevin equations (1), we employ an explicit first-order method as described in Ref. [48], where the boundary conditions can be accounted for by adding a term that ensures normal reflections at the boundary [49].However, to make the effect of the boundaries explicit in analytical calculations of the dissipation rate, we mostly employ smooth potentials that approximate the no-flux boundary conditions of the system.The potentials are where (a, m) ∈ Ω and the function V b with arguments a or m is given by The potential V b is a superposition of exponential functions which approximate hard, reflecting walls in the limit of small σ > 0. Non-zero values of σ ensure that the potential is smooth and bounded for r ∈ (0, r max ).To test if the finite range of the wall potentials affects the simulation results, we also perform simulations using "reflection operators", the details of which can be found in the Supplementary Material.The results of computations using the reflection operators are also presented in the Supplementary Material (FIG.S2) and agree with results that were obtained from computations that make use of the boundary potentials.

D. Analytical and Numerical Results
In this section, we study how the finite variable-range in a fully adaptive system leads to a non-zero adaptation error that increases with the magnitude of the stimulus.Concomitantly, the boundaries reduce dissipation, which results in a novel functional form of the cost-accuracy trade-off.
Since the system consists of fluctuating state variables, understanding its adaptive performance requires a consideration of the probability distribution of the variables.The probability distribution is affected by the presence of phase-space boundaries that determine the range of the system variables.We define the bulk of the phase space, denoted by B ⊆ Ω ⊂ Ω, as the largest subset of the phase space in which the boundary-induced changes of the probability distribution do not result in measurable changes in adaptation error.
The capacity of a system to adapt after a change in stimulus s is limited by the phase-space boundaries that determine the maximum and minimum values of system variables [50,51].In the present model, this claim can be proven by separating the adaptation error, Eq. ( 12), into three contributions, with which are obtained using the time-scale separation of adaptation and response.In Eq. ( 20), ϵ 1 and ϵ 2 are contributions due to the phase-space boundaries, while ϵ 3 accounts for the error in the bulk.For perfect adaptation, β = 1, the bulk contribution vanishes [38].Therefore, in a fully adaptive state, any error in the stationary value of the adapted state results from the limited range of the internal variables.Since in a perfectly adaptive state the adaptation error is only non-zero near the boundaries, the adaptation error for a large stimulus s depends on the maximum value that the methylation variable can assume.This maximum value is determined by m 0 .Results from numerical computations displayed in FIG.2a show that the dependence of the adaptation error on m 0 can be fitted by a Fermi-Dirac-like function as where c 1 , c 2 = c2 k s ,and c 3 are fit parameters and k = 10 4 is a constant scale factor.In FIG.2a, the values of m 0 are shown at which the analytical solution of the deterministic system predicts an error of 5% for different s.These values indicate the points where the boundary effects start to become significant for a given stimulus.Figure 2b illustrates that the adaptation error also depends on the setpoint a 0 .A larger setpoint results in an increased adaptation error.

Simulation Results: Dissipation Rate
If the stimulus magnitude s exceeds a threshold s 0 , the solution of Eq. ( 9) yields a fixed point m * > m 0 that lies outside the variable domain, m * ̸ ∈ Ω.Then, the FIG. 2. (a) Dependence of the adaptation error on the range of the control (methylation) variable m ∈ [0, m0], with k = 10 4 .For large m0, the error vanishes for a fully adaptive system since the state rarely encounters the boundary at m0. Simulation results can be fitted with an exponential form, see Eq. (22).Vertical lines show values of m0 at which the corresponding deterministic systems produce an error of 5%.(b) Adaptation error as a function of stimulus strength at fixed system size with m0 = 4. (c) The dissipation rate Ẇdiss of the adaptive system decreases with reduced system size if the state encounters the phase-space boundary.For small stimuli, s = 10 2 , the state remains in the center of the space and the dissipation is is independent of system size.(d, e) The function ∂mG(s, m)| m=⟨m⟩ss and the dissipation rate are roughly quadratic functions of the setpoint distance from the neutral state a0 = 0.5, which lies in between the active state (a = 1) and the inactive state (a = 0) of the system.(f ) Dissipation rate as a function of the stimulus strength s.Simulation results Ẇdiss are compared with an analytical approximation Ẇ a diss and an empirical formula W(s, λ).The three quantities agree for small stimuli s < s0.Otherwise, Ẇ a diss differs from the simulation results due to the breakdown of the analytical approximation at the boundaries.
steady-state expectation value ⟨m⟩ ss deviates from m * .To compute the dissipation rate Ẇ a diss , see Eq. ( 15), outside the bulk B, we replace m * with ⟨m⟩ ss .A comparison of m * and ⟨m⟩ ss as a function of the stimulus s can be found in the Supplementary Material, FIG.S3.
The interaction of state trajectories with phase-space boundaries reduces the system's energy dissipation rate in the stationary state.As shown in FIG.2c, for a small value of the input stimulus, s = 10 2 , the dissipation rate of the system remains constant for varying m 0 .However, for larger s the dissipation of the system starts to decrease when the probability distribution touches the boundary at m 0 .Furthermore, when s is fixed, changing the setpoint a 0 of the activity variable also affects the energy cost, as moving the distribution along the a−axis introduces effects of the boundaries at a = 0 or a = 1.The plot in FIG.2e shows that the dissipation rate decreases symmetrically as a 0 is moved away from the center at a 0 = 0.5.

Breakdown of Approximation for the Dissipation Rate and Empirical Formula
Comparison of the dissipation rate in simulations, Ẇdiss , with the approximate expression Ẇ a diss , Eq. ( 15), shows that the latter only holds if the setpoint is close to the center of the phase space, a 0 = 0.5, see FIG. 2e.Thus, the analytical approximation holds only if the effect of the boundaries is very small whereas Ẇdiss ≤ Ẇ a diss otherwise.
The behavior of the dissipation rate Ẇdiss can be further analyzed by defining a threshold value s 0 for the input stimulus.Above this threshold, the probability distribution is sufficiently altered by the boundary at m = m 0 to change the dissipation rate by more than some fixed percentage.For a 0 ∈ (0.3, 0.7), see FIG. 2e, and a large enough value of s such that m * ̸ ≈ 0, we conjecture the following expression for the dissipation rate  23) describes the numerical data well, while the approximation Ẇ a diss fails.

Cost-Accuracy Relation
What is the energetic cost of accurate adaptation?A well-known energy-speed-accuracy relation, proven in Ref. [38], asserts a monotonous increase of the dissipation rate with adaptation accuracy for varying chemical driving force.More precisely, the energy-speed-accuracy relationship states that the dissipation rate is proportional to the negative logarithm of the error.See the Supplementary Material (FIG.S4) for a recapitulation of this result in our simulations where the chemical driving is controlled by the value of β.Optimal accuracy is reached in the limit of perfect adaptation where the dissipation is maximal.
To understand how the finite range of internal variables affects the relationship between dissipation and accuracy, we vary here the magnitude of the external stimulus s.We focus on the limit of a fully adaptive model, i.e., β = 1.As shown in FIG.3a, the adaptation error obtained in simulations decreases linearly with a growing dissipation rate.In other words, the accuracy increases with the energetic cost also if the error solely results from a finite variable range.This relationship is thus qualitatively similar to the energy-speed-accuracy relation, but the functional forms and implications are different.
To buttress our numerical result for the cost-accuracy relation, we next perform analytical calculations of the dissipation rate.As illustrated in FIGs.2e and 2f, the validity of the approximate expression Ẇ a diss breaks down when significant adaptation errors occur close to phasespace boundaries.Thus, we derive a new equation for the dissipation rate by explicitly taking into account the boundaries of the variable domain.We consider the system of Langevin equations with any smooth respulsive boundary potential V as These equations are linearized around any point near the phase-space boundary (a c , m c ) of their deterministic counterparts.Using Eqs.(4, 5) and β = 1, we obtain where the higher order terms of the expansion are ignored and the leading terms in the Langevin equations drop out, as we define the fixed point (a c , m c ) to satisfy, The above two equations represent a modified version of Eq. ( 9).The boundary potentials ensure that the solution m c exists for all values of the input and that m c = m * for s < s 0 .The fixed point (a c , m c ) corresponds to the point of maximum likelihood -the point in the phase space where the probability density function reaches its maximum value, see FIGs.3b and 3c.The FPE corresponding to the linearized equations ( 25) is solved by a Gaussian ansatz.Substitution of the solution into the expression for the dissipation rate, Eq. ( 14), yields Here, ∂ 2 a V and ∂ 2 m V are both non-negative due to the choice of boundary potentials.The superscript "e" in Ẇ e diss indicates that this expression arises from a linear expansion of the governing equations of the system.
The approximation in Eq. ( 27) is only valid when either s < s 0 or s → ∞ because the state probability distribution can be approximated by a Gaussian distribution only for such values of s, see FIG. 3b.For other values of s, the distribution is skewed and an expansion around the point of maximum likelihood results in larger errors.
As illustrated in the Supplementary Material, simulations with the chosen exponential boundary potentials yield very similar results to simulations employing a reflection operator with very small timesteps.Thus, the results are expected to also hold for reflective boundaries.Also, a master-equation based chemotaxis receptor model with discrete molecule states yields qualitatively the same cost-accuracy relation for varying stimulus strength as shown for the Langevin-equation based model in this section, see Supplementary Material.

E. Discussion of Results for a Single Sensory System
Using a minimal Langevin-equation-based model for adaptation [38], we analyze how a limited range of internal variables affects adaptation accuracy and energetic cost.The model consists of a stochastic variable a called activity, a slow, stochastic control variable m, and an externally controlled stimulus s.In this work, the adaptation error ϵ is defined as the systematic, relative deviation of the activity a from its setpoint a 0 , see Eq. ( 12).The adaptation process can be visualized by following the motion of the system state (a, m) in the rectangular variable domain Ω with a ∈ [0, 1] and m ∈ [0, m 0 ].For the fully adaptive model considered here, the adaptation error solely results from the limited range of the internal variables a and m.To explain the onset of the error ϵ, consider the state to be initially at its setpoint with ⟨a⟩ = a 0 .After a sudden increase of the stimulus s, ⟨m⟩ slowly shifts to higher values to bring ⟨a⟩ back to its setpoint.However, if s becomes larger than a certain inputthreshold s 0 , such that the control variable approaches its upper limit m 0 , the reflective boundaries alter the probability distribution, see FIGs.3c and 1d.Then the control variable cannot increase sufficiently to counteract the effect of the external stimulus on the activity a and the perfect adaptability of the system is compromised.Unlike the error due to the limited range of the control variable, which results from the deterministic dynamics of the system, the error resulting from finite domain of a is a consequence of the finite, non-zero variance of the noise in the system.Also, as shown in FIG.2a, the adaptation error saturates for very large stimulus s to a maximum value lim s→∞ ϵ ≡ ϵ max < 1. ϵ max is independent of both the setpoint and the system size as long as m 0 is finite.The limit on the maximum error results from the skewness of the probability distribution.
The considered energetic cost of adaptation is the steady-state dissipation rate after a change of the input signal.Perpetual dissipation in the adaptive system results from the non-equilibrium forces that maintain the state (a, m) close to a fixed point in the presence of fluctuations.The fluctuations drive the state away from its fixed point and work is required to bring it back.Upon strong stimulation, the probability distribution for (a, m) moves to the vicinity of the boundaries.Intuitively, the effect of boundaries can be interpreted here as passive forces that also maintain the state at a fixed point.More formally, Eq. ( 10) ensures that the normal component of the fluxes at the boundary vanishes, thereby reducing the contribution of these fluxes to the integrand in Eq. (14).Consequently, the dissipation rate decreases at the boundary, see FIG. 2f.
A plot of adaptation error against dissipation rate as a function of stimulus strength yields a linear cost-accuracy relationship as a result of the limited range of internal variables.Thus, our result pertain to the behavior of a sensory system in a regime where its internal variables reach their physical limits.In that case, we see that the system is unable to adapt, but it also reduces the energy it uses for adaptation in a linear fashion.By contrast, variation of the chemical driving force yields a logarithmic cost-accuracy relationship, called the energy-speedaccuracy relation [38].Thus, the functional form of costaccuracy relations depends on the nature of the internal constraints imposed on the adaptation process.We emphasize, however, that the energy-speed-accuracy relation can directly imply a design trade-off for cells since the adaptive performance can be improved by investing more energy.This is not the case for the cost-accuracy relationship studied here, which depends on the range of internal variables and the external stimulus.In order to turn our results into a proper trade-off relationship, one would need to determine how the cost of enlarging the state space, i.e., to generate molecules with more methylation sites in the case of chemotaxis, improves chemotaxis.

III. INTERACTING SYSTEMS
We now shift the focus from single sensory units to pairs of interacting sensory units as illustrated in FIGs.4a.(i) and 4b.(i).The state of a sensory unit with index i ∈ {1, 2} is described by the variables (a i , m i ) that represent the activity and the control (methylation), respectively.The four-dimensional phase space is denoted as Ω

A. Overall Adaptation Error
Before specifying how the units interact, we introduce the adaptation accuracy of a combined system generically.For the coupled pair of systems, the overall activity a and overall error ϵ are defined as where ⟨a⟩ ss represents the stationary-state average of the overall activity and a 0 is the setpoint value, which is as- sumed be the same for both units.By employing the definition of overall activity, Eq. ( 28), and the sub-additivity of the norm, we can establish an upper bound on the overall adaptation error for the interacting system Here, ϵ 1 and ϵ 2 are the errors of the individual units, which are calculated with respect to the joint probability distribution of the interacting system.The upper bound (30) can be generalized to a system comprising an arbitrary finite number of sensory units.

B. Coupled Activities
The first interacting system we investigate involves coupling of the activities of two units.Coupling is as-sumed to be linear in the variables with symmetric coefficients.

Model Description
The system is described by a set of two Langevin equations for each unit, denoted with the indices i, j ∈ {1, 2}.For unit i and with neighboring unit j ̸ = i the equations are given by ȧi = Fai (a i , a j , m i , s) where η {mi,ai} is Gaussian white noise as in Eq. ( 1) and all η ... with subindices {a 1 , a 2 , m 1 , m 2 }, are statistically independent of each other.The generalized force Fai that appears in Eqs. ( 31) is defined as As the interaction between the two units is assumed to be symmetric, we have ω 12 = ω 21 .In Eq. ( 32) we set which represents the driving forces that depend on the input signal s through The time evolution of the joint probability distribution function P ≡ P (a 1 , m 1 , a 2 , m 2 , t) is governed by a Fokker-Planck equation.Conditions for detailed balance in absence of driving can be derived analogously to the non-interacting case, see Supplementary Material.These conditions are consistent with the assumption of one effective temperature for the two units with identical friction and noise strength coefficients.To control the interaction strength, we introduce a tuning parameter k such that where 0 ≤ k < ∞ and ω a1 = ω a2 .Furthermore, the detailed-balance conditions require that F mi is independent of a j for i ̸ = j and i, j ∈ {1, 2}.Integration of the detailed-balance conditions allows one to systematically construct an expression for F mi as where 0 ≤ β ≤ 1 and = 1 since the effective temperatures are assumed to be equal.We take β = 1 unless specified otherwise.
The boundary conditions for the system in Eq. ( 31) are similar to Eq. ( 10), but the probability flux is now a fourdimensional vector denoted by J = (J a1 , J a2 , J m1 , J m2 ) T .The boundary potential is re-defined accordingly.
The dissipation rate for the interacting system can be obtained from the dissipation rate integral of the system, analogous to Eq. (14).For system consisting of a pair of sensory units we obtain

Numerical Results for Coupled Activities
The set of Langevin equations (31) are solved numerically with the Euler-Maruyama method, as for the non- Figure 5a.(i)illustrates the relationship between the overall adaptation error and the dissipation rate for different stimulus magnitudes s.As for single sensory units, see FIG. 3a, the cost-accuracy relation resulting from the numerical solution of Eq. ( 31) is linear.For vanishing interaction strength, k = 0, the pair of units together have the same adaptation accuracy as an individual system, albeit with twice the dissipation rate.As k is increased, the slope of the linear relation becomes more negative, implying that interaction reduces the dissipation rate.Figure 5a.(iii)illustrates that a linear coupling of the activities reduces accuracy of the system.However, the overall adaptation error increases only slightly while the dissipation is strongly affected by coupling.

C. Coupled Methylation
After studying coupled activities, we next consider a coupling the control (methylation) variables.Again, the interaction is assumed to be linear and symmetric.

Model Description
In analogy to Eq. ( 31), the dynamics of a pair of sensory units labeled with the indices i, j ∈ {1, 2} and with coupled variables m 1 and m 2 is described by Langevin equations.For unit i and with neighboring unit j ̸ = i where the white noises η {ai,mi} are again pairwise independent of each other.The generalized forces F ai are the same as in Eq. (33).To guarantee detailed balance in equilibrium, the system must satisfy six constraints as shown in the Supplementary Material.These constraints are used to construct F mi , resulting in the same expression as given in Eq. (35).The coupling constants are scaled with the effective temperature T as where the parameter 0 ≤ k < ∞ is again used to tune the interaction strength.Lastly, as the interaction strength is symmetric and linear, the dissipation rate can be again computed using Eq. ( 36).

Numerical Results for Coupled Methylation
Results from simulations of the system described by Eq. ( 37) are shown in FIGs.4b(ii)-4b.(v)and 5b.The former displays the histogram of trajectories for different variable combinations, representing the probability distribution on different two-dimensional planes.The latter illustrates the cost-accuracy relation of the interacting system for varying interaction strengths.
Typically, the range of the control variable is much larger than the amplitude of the noise, m 0 ≫ ∆ m .Therefore, the probability distribution in the m 1 − m 2 plane appears highly localized, see FIG. 4b.(v).Coupling the control elements reduces the spread of this distribution on the m 1 − m 2 plane even more.The two slow variables are strongly correlated.However, the two fast activities are almost independent of each other and their joint distribution is only slightly elongated along the a 1 = a 2 axis, see FIG. 4b.(iv).Overall, the probability distribution on the a 1 − m 1 and a 2 − m 2 planes are very similar to the distributions for a non-interacting system, FIGs.(4b(ii), 4b.(iii)).Therefore, the distributions of the coupled units are distorted after application of a stimulus almost the same way as for a non-interacting system.Adaptation errors of the fully adaptive, interacting system are consequently almost the same as for a single sensory system.
A coupling of the control variables reduces the overall dissipation rate both for trajectories that stay in the phase space interior and for trajectories that interact with the boundaries.For a given coupling strength k, the dissipation rate is slightly higher if the m i are coupled, compared to coupled activities a i .Figure 5b.(i)shows a linear relationship between the adaptation error and the dissipation rate required for adaptation, similar to FIG. 5a.(i).Overall, coupling of the control variables reduces the dissipation rate and the adaptation error appears to be independent of interaction strength for this coupling, see FIG. 5b.(iii).

D. Discussion of Interacting Systems
To study the adaptive performance for an idealized system consisting of two sensory units, we define an overall adaptation error in terms of the overall activity, see Eqs. (29,28).These definitions result in an upper bound for the overall adaptation error of the system in terms of the individual errors of the interacting sensory units.If one measures the error of individual adaptive units, the overall error must be smaller than or equal to the mean of the individual errors.This inequality is generic as it is robust to the type of interaction between the units and is true even if the systems are not identical, provided they have the same setpoint a 0 .
The first interacting system that we consider is one in which the activities a i are coupled, see Eq (31).For this system, we find that increasing the interaction strength leads to a reduction of the dissipation rate and a slightly worse adaptation accuracy.This change of the system behavior is due to a stretch of the probability distributions for the (a i , m i ) along the m i −axes.After an increase of the stimulus s, the probability distributions shift to higher values of the control variables m i .During this shift, the distributions interact earlier with the boundaries of the phase space if they are stretched.Thus, boundary conditions of the variables influence the behavior of activity-coupled units at smaller stimuli s than for the corresponding non-interacting unit.Instead of coupling the fast activities, we next couple the slow control variables, see Eq. (37).Similarly to the previous case, the linear, symmetric coupling of the control variables reduces the dissipation rate.However, the state probability distributions for the individual units in the system closely resemble the distribution for a single, uncoupled unit.Therefore, the adaptation accuracy is not compromised by this coupling.
For both types of interacting systems, the finite range of internal variables is irrelevant for small stimulus values s if the states (a i , m i ) rarely assume their extreme values.In absence of such boundary effects, an increase of the coupling constant reduces the overall dissipation from originally twice the dissipation rate of a single unit down to almost the dissipation rate of only one unit.Thus, a linear interaction of sensory units in the fully adaptive limit reduces the number of wasteful cycles and effectively lowers the energetic cost of adaptation.

IV. CONCLUSION
From bacterial chemosensors to yeast osmotic pressure sensors and mammalian light sensors, a wide variety of systems in biology rely on biochemical feedback networks to enable sensory adaptation [2,52,53].In the case of E. coli chemotaxis, transmembrane receptors sense and adapt to the concentration of extracellular ligands.This process is characterized by a rapid response in sensor activity following a change in ligand concentration and a slowly varying methylation that counteracts the effect of the stimulus and returns the activity to a set point [32].
We employ an idealized model of bacterial adaptation to visualize how the working range affects the accuracy and energetics of adaptation.First, we study the threenode model for adaptation of a single sensory unit.Large stimuli drive the state of the modelled system to its physical limits, represented by boundary conditions of the phase space of internal variables.We find that the adaptation error is a linear function of the dissipation rate.As the strength of a stimulus is increased, the error increases while the dissipation rate decreases.The results from simulations are supported with analytical approximations.Next, we investigate the the same relationship for interacting systems.For pairs of sensory units, the overall adaptation error is determined by the average activity of the two units.An interaction between the units is modelled in two ways.Either by coupling the output activities of the two units or by coupling their control variables representing the methylation state.We find that both types of interactions reduce the dissipation rate and therefore affect the cost-accuracy relationship.Furthermore, the simulation results suggest that a coupling of the methylation variables may be more advantageous than a coupling of the activities because, in this case, the interaction reduces the energetic cost cost of powering multiple units without compromising the accuracy of adaptation.Overall, a cooperative operation of many sensors amplifies sensitivity [21] and our models predict that the amplification can be energy efficient since the dissipation rate per sensor is reduced by coupling.
Looking ahead, one can consider to study the energetics of larger systems of coupled units, as well as how they behave in the regime of imperfect adaptation.Clearly, linear receptor arrays or randomly clustered receptors with next-neighbor couplings may exhibit more complex collective energetics.From a more formal point of view, it would also be interesting to study how different levels of noise in the different components of the system would affect the adaptive performance while at the same time driving energy flows.Finally, further theoretical and experimental work is needed to assess the biological relevance of a trade-off between the benefit of a large working range and the energetic cost of generating sufficiently complex sensory systems.The analysis of the scenarios involved is far from straightforward, but this work represents a step towards understanding the energetic consequences of a limited working range in adaptation.

FIG. 1 .
FIG. 1.(a) A minimal three-node network with a negative feedback loop.(b) Following a change of the stimulus s, the instantaneous response of the activity a is cancelled by a gradual change in the control variable m.The adaptation error is defined as the deviation of the expectation value of a from its setpoint a0 in the stationary limit.(c) Phase portrait of the dynamics in the absence of noise.For β = 1, the fixed point of the deterministic dynamics (a0, m * ) is stable and represents the adaptive state of the system.(d) Steady-state probability distribution of the system variables in the presence of noise.The deterministic fixed point lies at the center of the probability distribution, that is, (⟨a⟩ss, ⟨m⟩ss) = (a0, m * ).(e) For β = 0, the chemical driving vanishes and the fixed point (a0, m * ) is an unstable saddle node.(f ) In the presence of noise, the state of the non-adaptive system is driven to the corners of the phase space.Simulation parameters: s = 20, ωa = 50, ωm = 5, T = 0.01.For (f ), we chose s = 3 to illustrate that the state localizes in both corners of the phase space.
) where K 0 ≡ Ẇ a diss s=s0 is the dissipation coefficient and ϵ 0 ≡ min s∈(0,∞) Ẇdiss .Results from the empirical expression W diss and Ẇdiss are shown in FIG.2f.The expression in Eq. (

FIG. 3 .
FIG. 3. (a) Cost-accuracy relation for varying stimulus strength.The adaptation error ϵ is linearly related to the the dissipation rate Ẇdiss at different values of s.Analytical results Ẇ e diss estimating the dissipation rate at weak and strong stimuli, s < s0 and s ≫ s0, agree with the simulation data.The steepest decent approximation used for Ẇ a diss breaks down for strong stimuli.(b) Steady-state probability distribution after application of a weak stimulus, s = 10 3 < s0.The distribution resembles a Gaussian and the averages co-localize with the solution of Eq. (26), (a c , m c ). (c) State probability distribution after application of a strong stimulus, s = 5 • 10 4 > s0.The distribution is skewed and the averages do not co-localize with the maximum of the distribution (a c , m c ). Simulation parameters: V0 = 1/2, ωm = 5, ωa = 50, σ = 0.01, T = 0.01.

FIG. 4 .
FIG. 4. (a.(i), b.(i)) Composite systems consisting of two sensory units, where each realizes an adaptive feedback loop.The units are linearly coupled either via their activity or their control (methylation) variables.The overall activity a of the system is the average of the individual activities a1 and a2.(a.(ii) -a(v)) Steady-state probability distribution P visualized on two-dimensional planes in the phase space for coupled activities.A large coupling constant ω12 > ωa {1,2} leads to a strong correlation of the rapidly changing activities a1 and a2.The control variables m1 and m2 are also strongly correlated.As a result, the distributions on the ai − mi planes for the individual units are also altered.(b.(ii) -b.(v))State probability distribution P visualized on two-dimensional planes in the phase space for coupled control variables.Due to the slow dynamics of the control variables in the individual units determined by ωa {1,2} ≫ ωm {1,2} , the variances of the m {1,2} are smaller than the variances of the a {1,2} .As a result, even if a large coupling constant enforces a strong correlation of m1 and m2, the activities are only weakly correlated as seen by the almost radially symmetric distribution in the a1 − a2 plane.Consequently, the probability distributions for the individual units shown in (b.(ii) -b.(iii)) have largely the same shape as in absence of coupling, see FIG. 1d.Simulation parameters: k = 2000, s = 100, ωa {1,2} = 50, ωm {1,2} = 5, T = 0.01, m0 = 4.
interacting system.Results are shown in FIGs.4a.(ii)-4a.(v)and FIGs.5a.The probability distribution function P has four dimensions and plotting the distribution requires a selection of two-dimensional planes in the phase space as shown in FIGs.4a.(ii)-4a.(v).Figures 4a.(ii) and 4a.(iii) display distributions for two sensory units with strong coupling, i.e., k ≫ 1.Comparison with results for a single unit, see FIG. 1d, shows that the coupling distorts the probability distributions since it introduces a correlation between the two activity variables.Consequently, the probability distribution on the a 1 − a 2 plane aligns along the line a 1 = a 2 and the distributions in the two a − m planes also exhibit a stretched variance along the m−axis.

FIG. 5 .
FIG. 5. Cost-accuracy relation for composite systems consisting of pairs of linearly coupled sensor units.(a.(i)) Adaptation error vs. dissipation rate for varying stimulus magnitude s in systems with linearly coupled activities.As for a single, noninteracting unit, the cost-accuracy relation is linear.The parameter k controls the coupling strength.When k = 0, the units do not interact and the overall dissipation rate is twice the dissipation rate of a single unit.As the coupling strength k increases, less work is required to achieve the same accuracy as without coupling.(a.(ii), a.(iii)) Dissipation rate and overall adaptation error as a function of the stimulus magnitude for different coupling strengths k.With increasing coupling strengths, the dissipation rate decreases while the adaptation error increases slightly when activities are coupled.(b.(i)) Adaptation error vs. dissipation rate for varying stimulus strength s in systems with linearly coupled control (methylation) variables.For strong coupling, k ≥ 20k0, the overall dissipation rate nearly equals that of a single non-interacting system.(b.(ii), b.(iii)) Dissipation rate and overall adaptation error as a function of the stimulus magnitude for different coupling strengths k.With increasing coupling strength, the dissipation rate decreases while the adaptation error remains unchanged.Linear coupling of the slow control variables thus reduces the energetic cost of adaptation per unit without compromising the adaptive performance.(c.(i), c(ii)) The adaptation error is consistently lower if control variables mi are coupled, rather than activity variables ai, for otherwise identical systems.Dots are simulation results, dashed lines show data trend.Simulation parameters: k0 = 100, ωa {1,2} = 50, ωm {1,2} = 5, T = 0.01, σ = 0.01, V0 = 1/2.
the equations are given by ȧi