Complete integrability of information processing by biochemical reactions

Statistical mechanics provides an effective framework to investigate information processing in biochemical reactions. Within such framework far-reaching analogies are established among (anti-) cooperative collective behaviors in chemical kinetics, (anti-)ferromagnetic spin models in statistical mechanics and operational amplifiers/flip-flops in cybernetics. The underlying modeling – based on spin systems – has been proved to be accurate for a wide class of systems matching classical (e.g. Michaelis–Menten, Hill, Adair) scenarios in the infinite-size approximation. However, the current research in biochemical information processing has been focusing on systems involving a relatively small number of units, where this approximation is no longer valid. Here we show that the whole statistical mechanical description of reaction kinetics can be re-formulated via a mechanical analogy – based on completely integrable hydrodynamic-type systems of PDEs – which provides explicit finite-size solutions, matching recently investigated phenomena (e.g. noise-induced cooperativity, stochastic bi-stability, quorum sensing). The resulting picture, successfully tested against a broad spectrum of data, constitutes a neat rationale for a numerically effective and theoretically consistent description of collective behaviors in biochemical reactions.

Restricting to steady states, an alternative approach relies on statistical mechanics, as suggested by C.J. Thompson in his seminal work [26].Indeed, statistical mechanics turns out to be particularly effective for the description of universal behaviors of a wide range of biological systems (from the extra-cellular level of neural [27][28][29] or immune networks [12,30,31], to the intracellular level of gene regulatory and protein networks [30,32,33]).Moreover, as recently observed in the case of large systems [34,35], the statistical mechanical approach plays the role of a general stochastic framework that naturally highlights the structural and conceptual analogies between response functions in biochemical reaction kinetics and transfer functions in cybernetics (see Fig.s 1 and 2), thus tacitely working as a translator between these two worlds, that is crucial to show how information is handled by these biochemical systems.where the response is the output voltage while the stimulus is the input voltage.(b) The sigmoidal shape of the self-consistency of a ferromagnet, where the response is the magnetization and the stimulus is the external magnetic field.(c) The sigmoidal response of the activation function of a neuron where the output is the action potential voltage while input is conveyed by all the afferent electric synaptic currents.(d) The sigmoidal shape of the saturation curve of a cooperative chemical reaction, where the output is the fraction of bound sites, while the input is the (log-)concentration of the ligand.
However, the current research in biochemical information processing has been recently attributing particular importance to systems involving a relatively small number of units and this implies that the standard statistical mechanical picture, given in the thermodynamic limit (where the role of intrinsic noise can be suppressed [36]), is not accurate.One of the goals of the present work is to extend the theoretical framework developed in [34,35] to allow for the description of systems of finite sizes.
Before proceeding it is worth stressing that the statistical mechanical description of chemical kinetics pursued here is based on the canonical ensemble and it is accomplished at the mean-field level.This is therefore clearly different from the statistical mechanical description based on the gran-canonical ensemble introduced in the 70's [37].The advantage of the present formulation is that the use of a canonical framework allows us to establish structural bridges between the phenomenology of these bio-chemical systems on one side and the well-consolidated theory of information processing systems (i.e.cybernetics) on the other side, since the latter (inflected for instance in terms of neural networks and learning machines) is mainly developed within the canonical formalism [28].Along the same line we choose to adopt the mean-field perspective: admittedly, this implies a cost (as we give up a detailed description of the true architecture of the system under consideration and we deal with effective parameters to be properly renormalized), yet the reward lies in the possibility to directly compare the emerging response functions with transfer functions in cybernetics and therefore to understand how information is processed through a given reaction.Further, trough direct calibration of the re-normalized key parameters over a small subset of data, the whole theory become of immediate experimental applicability.
Let us now present in more detail the statistical-mechanical framework adopted in this work and the results we obtain through this path.We consider large macromolecules (e.g., polymers, proteins, etc.), whose binding sites are dichotomic variables (Boolean logical operators to be thought of as Ising spins) that can be either empty or occupied by a ligand (e.g., a substrate): the log-concentration of free ligands plays the role of the external field in spin systems.The binding sites are not necessarily independent: in the special case where they are independent the emerging kinetics follows the so-called Michaelis-Menten law, and this is naturally captured by the lack of interaction between spins within the statistical mechanical route of formalization.On the contrary, if there is interaction, the kinetics can be cooperative (for positive values of the interaction), or anti-cooperative (for negative values of the interaction), and these behaviors are typically coupled to various names, e.g., Hill, Koshland, and Adair reactions, etc. (see e.g.[37]).More precisely, binding sites with identical structure/function (e.g., the four oxygen-capturing arms of the hemoglobin) are modeled as a unique spin system, where each spin feels the external field (the oxygen log-concentration in this example) and interacts imitatively with other spins; conversely, binding sites with different ) obeys cooperative kinetics, the second (calmodulin dependent protein kinase 2) ultra-sensitive kinetics, while the last (synaptic glutamate receptor) fulfills anti-cooperative kinetics.In the second row the saturation curves for these proteins are shown: α is the concentration of the free ligand needed for their binding and Y is the fraction of occupied binding sites.Symbols with the relative error-bars stand for real data taken from [20,38,39] and lines are best fits performed through the analytical expression in Eq. ( 17), where the best-fit coefficient J corresponds to an ergodic ferromagnet (left), a low-temperature ferromagnet (center) and an anti-ferromagnet (right).In the third row three circuits are shown.From left to right, an operational-amplifier, an analogto-digital converter and a flip-flop, while in the fourth row their transfer functions are presented to highlight the behavioral analogy with the proteins.The latter shows the response of the system (output voltage) as a function of the input of the system (input voltage).Following columns instead of rows, the first column is due to cooperative behavior, the second to ultrasensitive behavior and the last to anticooperative behavior.The authors are grateful to Nature Publishing Group for the permission to reproduce this figure from [34].
structure/function are modeled as a spin system fragmented into more parties, where spins belonging to different parties compete to bind the ligand by anti-imitative interactions (e.g., for the two insulin-capture α-subunits in the insulin receptor, the system describing the receptor is split into two main parties -i.e.arms-that anti-cooperate).This scheme naturally leads to a statistical mechanical scenario consisting in multi-partite spin systems (on which -in turn-there has also been recent interest in the Statistical Mechanical community [40][41][42][43][44][45]), with (positive or null) intra-party interactions and (positive or negative or null) inter-party interactions; we will focus on the two-party case, as schematized in Fig. 3.The behavior of this system in the presence of an external magnetic field h (i.e., the input) can be captured in terms of the average magnetization m (i.e., the output), which keeps track of the underlying collective features among spins.From a kinetics perspective, the function m(h) recovers the saturation function of a system of binding sites as the concentration of free ligands is tuned.The isotherm of the magnetization and the saturation function are just two different inflections of the same transfer function, namely of the same input-output relation.In this paper the exact, explicit, solution of these magnetic systems (and of the reaction kinetics they code for too) is given also at finite volumes.In particular, we derive a set of linear multidimensional partial differential equations for the partition function of the bipartite spin model which are valid for any (either finite or infinite) number of spins N .The solution for the model is specified via a suitable initial datum.We show that in the thermodynamic limit N → ∞ the free energy related to these systems fulfils a set of Hamilton-Jacobi type equations that is completely integrable by the characteristics method and we provide the explicit solutions in terms of partial magnetizations, which in turn satisfy a set of two coupled equations of state.Via this route, the way these reactions handle information processing becomes transparent; such equivalence is summarized in Fig. 2   The cooperative case is characterized by a Hill coefficient nH larger than 1, corresponding to a positive coupling J between spins.The anti-cooperative case is characterized by a Hill coefficient nH smaller than 1, corresponding to a negative coupling J between spins.The input of the system is provided by the concentration of free ligand α, corresponding to the magnetic field h (more precisely, h = 1/2 log(α/α0), where α0 is the half-saturation concentration, in such a way that a negative field is equivalent to a concentration smaller than α0 and vice versa).The output is provided by the saturation Y , corresponding to the magnetization m (more precisely, m = 2Y − 1, in such a way that the mean performed over the magnetizations of the two parties is equivalent to the sum of the saturations pertaining to the two binding sites minus 1).The plots in the bottom represent the evolution of the output as the input is varied.In particular, the curves shown here are obtained as the solution of Eq. ( 17) for J = 0.5 and J = −1.2,respectively.
large systems: we show that the key parameters of the theory match -always with remarkable accuracy-the standard empirical indicators (e.g., the coupling constant used in the Hamiltonian spin-like representation of the reaction recovers the Hill coefficient of standard kinetics literature).Then, we proceed by considering small systems (i.e.away from the thermodynamic limit) and we show that the solution at finite sizes can still be obtained explicitly (by separation of variables once an Ansatz on the initial state is given): the framework obtained in this way is used to explore and address several phenomena recently highlighted in the experimental literature on small system's kinetics (e.g., we recover that systems devoiding of cooperativity can still display a cooperative-like behavior and, possibly, bistability phenomena due to stochastic effects, as for instance discussed in [46]).
The paper is organized as follows.First, we provide a brief survey of the statistical mechanical formulation of reaction kinetics as proposed by Thompson in [26] and later developed in [34,35].Then, we derive differential identities for the partition function and the free energy related to these models and construct their solutions both in the thermodynamic limit and for small systems.Next, we discuss implications of our theory such as noise-induced cooperativity and bistability, signal amplification and noise suppression.Finally, we summarize and comment on our results and discuss further outlooks.

Standard chemical kinetics: statistical mechanical formalization and cybernetic interpretation
Hereafter we review main concepts of chemical kinetics, statistical mechanics and cybernetics, which we will be needed in the following; we refer to classical textbooks for a more extensive treatment (see e.g., [36,[47][48][49][50]).
Chemical-kinetics framework.In many macromolecules (i.e.polymers, proteins, etc.) ligands bind in a nonindependent way.In particular, if, upon a ligand binding, the probability of further binding (by other ligands) is enhanced, as for example in the paradigmatic case of hemoglobin [26], the system is said to exhibit positive cooperativity; viceversa, cooperativity is negative where further binding of ligands is inhibited [51] as in the case of some insulin receptors [52] (see Fig. 3 for a schematic representation).Fundamental mechanisms underlying cooperativity depend, in general, on the microscopic details of the system under consideration.For instance, in the case of polymers, if two neighbor docking sites bind charged ions, the electrostatic attraction/repulsion may be responsible for positive/negative cooperativity.
In complete generality, let us consider a model where several identical hosting (macro)-molecules P can bind overall N identical small molecules S (whose concentration is denoted by [S] ≡ α) on their structure; calling P j the complex of a molecule P with j ∈ [0, N ] molecules attached, at the chemical equilibrium we have For the sake of simplicity, let us consider the case j = 1.The time evolution of the concentration [P 0 ] of the unbound protein P 0 is governed by the equation where K +1 , K −1 are, respectively, the forward and backward rate constants for the state j = 1, and their ratio defines the association constant −1 .In the steady state d[P 0 ]/dt = 0, we have For generic j > 0 we have Therefore, in general, one can write , and, by extension, [P j ] = [P 0 ]( j i=1 K (i) )α j .In many practical situations, a direct measure of [P j ] is not feasible.A convenient experimental observable is then given by the average number S of ligands that, at the equilibrium, are bound to the macromolecule(s) and this is given by The last expression in the previous equation is the well-known Adair's equation [36], obtained by iterating the relation (2).In this work (as standard) we will focus on the normalized version of S, called saturation function Y and defined as Unless otherwise stated, Y represents the expected fraction of occupied sites in the whole set of molecules P under investigation, namely the global system's response (i.e. the output signal) to the stimulus provided by the ligand's concentration α (i.e. the input signal).In a non-cooperative system, one expects independent and identical binding with microscopic association constant K, and one can write K (j) = (N − j + 1) K/j.In this case, Adair's equation (3) for S and Y reads as The latter expression is the well-known Michaelis-Menten equation [36].
Clearly, the kinetics becomes far less trivial as soon as interactions among binding sites occur.For the sake of simplicity, if we consider the limiting case where intermediate steps can be neglected, i.e.
we get that More generally, accounting also for some degree of sequentiality [36], one obtains the well-known Hill's equation where n H , referred to as Hill coefficient, represents the effective number of interacting binding sites.We note that, by posing n H = 1, equation ( 9) recovers the Michaelis-Menten law (6).In fact, in the non-cooperative case, the number of binding sites effectively interacting is just one (i.e., there are no true interactions).For n H > 1 the kinetics is said to be cooperative, while for n H < 1 it is said anti-cooperative.If n H 1 the kinetics is said to be ultra-sensitive.Given a set of experimental measures of the saturation function Y (α), n H is estimated as the slope of log[Y /(1 − Y )] versus α, measured at half-saturation (namely, when Y = 1/2).
The statistical mechanics formalization.Statistical mechanics can provide a convenient language to describe biochemical kinetics and frame it into a (noisy) logical scaffold (see [34,35] for details).In the present section we briefly review how main concept of mean-field cooperative statistical mechanics (i.e.ferromagnetism) apply to the scenario of cooperative reaction kinetics.We start from the microscopic model described above, consisting of an ensemble of identical macromolecules, carrying overall N binding sites labeled by i = 1, 2, . . . .Macromolecules are in a solution with smaller molecules (the substrate) and each binding site can accommodate only one molecule of the substrate.Notice that, following the mean-field approach, here we do not distinguish binding sites belonging to the same molecule, but we are treating the whole set of binding sites, pertaining to the whole set of molecules, at once.
Let us associate to each binding site an Ising spin σ where σ i = +1 if the i th site is occupied and σ i = −1 if it is empty.A particular configuration of the system is specified by the set {σ}, through which all the classical observables can be introduced.For instance, the magnetization is defined as m({σ}) = ( N i=1 σ i )/N and is directly related to the occupation number S of binding sites by in such a way that the saturation function (in terms of the magnetization) reads as For simplicity, we first consider a system which does not exhibit collective behavior, that is, no interaction between binding sites is present and only the interaction between the binding sites and the ligand occours.The external field in statistical mechanics, namely a scalar parameter h, is introduced as the log-concentration log(α) of free-ligand molecules [26].As there are no spin-spin interactions, this system is mapped into a system of free spins σ thermalyzing in the presence of a magnetic field h and whose energy function (or Hamiltonian) is given by Note that h is assumed to be independent of the site index, meaning that binding sites interact uniformly and independently with the ligands; this also means that the time-scale for diffusion is fast with respect to the time scale for thermalization.
Once the Hamiltonian representing the model under investigation is defined, we can apply the standard statistical mechanics machinery, namely the Maxwell-Boltzmann probability distribution P ({σ}) associated to a generic system state {σ} where Z is the normalization, called partition function and β tunes the level of noise in the system, in such a way that for β → 0 the probability measure becomes flat and every state assumes the same chance to happen, while for β → ∞ the system collapses in configurations corresponding to the minima of the Hamiltonian.Throughout the paper, given a generic function of the spins f (σ), the bracket averages f (σ) represent the averages over the Maxwell-Boltzmann probability distribution.
An explicit evaluation of the partition function Z allows for a direct measurement of the free energy F related to the model encoded by H, that is, F = (1/N ) ln Z (notice that, in our derivation, just for mathematical convenience we work with the negative free energy F = −(1/N ) ln Z; hence, max(F ) just corresponds to min( F )).In general, F depends on the system configuration {σ} and on the system parameters, namely the external field, the level of noise, and the possible coupling between spins.Once these parameters are set, the maximization of F provides the related equilibrium state(s).This is because the free energy F is nothing but F = −βU + S, where U = H is the internal energy of the system (i.e. the Maxwell-Boltzmann average of the Hamiltonian at given noise β) and S = − ln P ({σ}) is the entropy, thus, maximizing F implies looking for the minimum of U and, simultaneously, the maximum of S (at given β).
In particular, for the system described by the Hamiltonian ( 12), the energy is minimized by those configurations where spins are aligned with the external field.Thus, in the absence of any source of noise, the system relaxes to a state where m = 1 (i.e., σ i = +1, ∀i) if h > 0, or to a state where m = −1 (i.e., σ i = −1, ∀i) if h < 0. In the presence of noise the equilibrium state requires, beyond energy minimization, also entropy maximization: setting for simplicity β = 1 (without loss of generality since this is equivalent to rescale h → hβ) one can find that the equilibrium state corresponds to m = tanh(h) (see e.g., [48]).Actually, one could reach the same result without relying on statistical mechanics, as explained hereafter.First, we need to relate h to the ligand concentration α.As mentioned above, the energy is minimized by those configurations where spins are aligned with the external field: if h > 0 molecules tend to bind to diminish energy, while if h < 0 bound molecules tend to leave occupied sites; this suggests that h plays the role of the chemical potential for the binding of the ligand molecules on the docking sites.Moreover, as observed in [26,34,35], the chemical potential can be expressed as the logarithm of the concentration of the substrate, upon proper normalization, namely where α 0 stands for the value of ligand concentration such that binding sites have the same probability of being occupied or unoccupied.
Crucially, under the mean-field approximation (and this is the only case, apart from linear chain models which, still, do not exhibit phase transitions and are of moderate interest), the probability P ({σ}) of the configuration {σ} can be factorized as P ({σ}) = N i=1 P (σ i ).One can therefore focus on the single-spin probability and, following [26], state that P (+1) is proportional to the concentration α and that P (−1) is proportional to the inverse of the concentration α, that is, P (+1) = Ce +h and P (−1) = Ce −h , where C is fixed in such a way that P (+1) + P (−1) = 1, i.e.C = 1/(2 cosh(h)).Then, the expression Exploiting this result the average saturation function (see Eq. ( 11)) reads as Using Eq. ( 13), and recalling that tanh , it is immediate to check that Eq. ( 14) coincides with the Michaelis-Menten equation ( 6) with the particular choice K = 1/α 0 .This is perfectly consistent with the underlying assumption of independent binding sites (i.e.no couplings among spins) tacitly made when we defined the system under study, via the Hamiltonian (12).
Let us now generalize this scenario by introducing the simplest possible two-body interaction given by the Hamiltonian The Hamiltonian (15) represents the well-known Curie-Weiss theory of ferromagnetism [48].The first sum in the right hand side of ( 15) accounts for all possible N (N − 1)/2 interacting pairs of spins and the coupling is homogeneous as the constant J is the same for all pairs.If J > 0 the configurations where spins are aligned are more favored thus this choice naturally leads to a theory for cooperative kinetics, while J < 0 favors the configurations where spins compete and are misaligned thus working for the anti-cooperative case.Clearly, models with different values of J represent different chemical systems.
As well known, the condition of minimization for the Curie-Weiss free energy in the thermodynamic limit N → ∞, yields the so-called self-consistency equation (see e.g., [48]) Recalling the mapping (11), and reabsorbing β by the rescaling βJ → J and βh → h, the previous equation translates into the reaction kinetics vocabulary as Eq. ( 16) implies that, at low noise levels, the Curie-Weiss model exhibits an abrupt change in the magnetization as a function of h.More precisely, a second order phase transition occurs at h = 0 (i.e., α = α 0 ) so that m, and then Y , is continuous while its derivative diverges as a function of J at the critical value J c = 1.The phase transition is first order if, at fixed J > J c , m (and then Y ), is viewed as a function of h.In this case, for J > J c , a jump occurs at h = 0 and the reaction is ultra-sensitive (its transfer function mirrors that of an ON/OFF switch).
As a robustness check, we verify that, in the absence of interaction, the model is consistent with the classical Michaelis-Menten kinetics.Indeed, Eq. ( 17) can be written as and, setting J = 0, the equation above clearly reduces to Eq. ( 6) with K = α −1 0 .In this case, no phase transitions occur and the model turns out to be not suitable to describe complex biochemical reactions.
As anticipated, in order to provide a quantitative information about how a system actually departs from the simplest Michaelis-Menten framework, the so-called Hill coefficient is introduced, defined as the slope of log[Y /(1 − Y )] versus α, which can be recast as where the last expression was obtained using Eq.(17).It is straightforward to check that the Hill coefficient for the Michaelis-Menten model, obtained for J = 0, is n H = 1 as expected.On the other hand, a positive value of J implies a positive cooperativity and the closer J is to the critical value J c = 1 (from below) and the stronger the cooperativity exhibited by the system; values of J larger than 1 correspond to discontinuous saturations (i.e.systems showing ultra-sensitive responses).We stress that the relation (19) provides a crucial link between the experimentally accessible quantity n H and the theoretical parameter J, namely the coupling strength of the underlying mean-field spin description.
We finally observe that small-coupling expansions of the right hand side of ( 17), or, equivalently, of (18), lead to suitable polynomial approximations for the saturation function to be possibly compared with formulas obtained through classical (usually ODE-based) phenomenological approaches in chemical kinetics [37].For example, the first order expansion of the expression (18) around J = 0 gives which is equivalent to Adair's equation (3) K (2) .The cybernetic interpretation.The cybernetic interpretation of chemical kinetics, extensively treated in [34,35], is based on the behavioral analogy between the saturation curves (or binding isotherms) in chemical kinetics, the self-consistencies (i.e.state-equations) in statistical mechanics and the transfer functions in electronics, see Fig. 2. Similarly to saturation curves, self-consistency functions and transfer functions are nothing but relations between an input (the ligand concentration α in chemical kinetics, the magnetic field h in statistical mechanics and the input voltage V in in electronics) and an output (the saturation function Y in chemical kinetics, the expected magnetization m in statistical mechanics and the output voltage V out in electronics).
As an illustrative example, let us consider the paradigmatic operational amplifier.In a regime of small input voltage V in , the output voltage V out , is described by the following transfer function where G = (1+R 2 ) is referred to as gain of the amplifier and R 2 is the feed-back resistor (allowing for true amplification as, if R 2 = 0, then G = 1), as shown in Fig. 2 (left column).Similarities between this response function and the saturation curve in chemical kinetics can be clarified, using the statistical mechanics formalism, by comparing the mathematical structure of the response of these systems, namely comparing V out in (21) with the average magnetization m of a ferromagnet in the linear regime [where x = (J m + h) is small such that tanh(x) ∼ x]: Observing that, for J < J c ≡ 1, the Hill coefficient ( 19) can be approximated as n H = 1/(1 − J) ∼ (1 + J), we can write The conceptual term-to-term identification between inputs and outputs of the above transfer functions suggests that the Hill coefficient can be interpreted as the gain factor for the reaction (and we can already see why cooperativity is needed to amplify biochemical signals, as we require J > 0).Let us also note that Hill coefficients are typically of order ∼ 10 or less, and this contributes to explain why amplification of biochemical circuits is rather low [53] if compared to its electronic counterpart, where the gain can be of several orders of magnitude [49].We emphasize that this behavioral analogy between biochemical kinetics, statistical mechanics and electronics goes far beyond the linear regime exploited above because all the related response functions (i.e.binding isotherms, selfconsistencies and transfer functions), far away from the linear regime, display intrinsic saturation effects (see Fig. 1. right panel): roughly speaking, if all the macromolecular binding sites are occupied, no matter how further we increase the ligand concentration, the system will not vary its response that is already maximal.The same applies to spin systems as, once all the spins are aligned with the magnetic field, a further increase in the latter can not produce any shift in the system's response and the same holds for operational amplifiers too as, once the output voltage reaches the collector tension, any further increase in the input voltage will no longer produce a variation in the response.Moreover, apart from the operational amplifier, also other devices naturally fit the picture provided for describing biochemical reactions in different regimes.For example, if J > J c (i.e., n H 1 and G 1) the expected magnetization m(h) as well as the saturation function Y (α) develop a discontinuity at h = 0 and α = α 0 respectively (and we refer to ultra-sensitive kinetics in biochemistry and first-order phase transitions in statistical mechanics).The corresponding limit for the operational amplifier is the analog-to-digital converter (ADCs) at V in = 0 (see Fig. 2, center column).Another example is given by the simplest bistable flip-flop, constituted by two saturable operational amplifiers reciprocally inhibiting, in such a way that the output of one of the two amplifiers is used as the inverted input of the other amplifier: this is the simplest 1-bit memory device as it is possible to assign a logical 0 (or 1) to one state and the other logical 1 (or 0) to the other state.As the two operational amplifiers (i.e. the two parties) interact inhibiting reciprocally, the translation of this circuit into a chemical circuit will require anti-cooperativity among the two parties, thus will be naturally accounted when anti-cooperativity is at work (see Fig. 2, right column).When dealing with more complex reactions, such as reactions involving several components and several steps, our approach allows to recognize the various basic 'devices' at work and to build the whole circuitry in a cascade fashion so to figure out how the equivalent 'electronic circuit' effectively processes information as the reaction goes by.In particular, a suitable combination of the elementary bricks described above leads to the construction of (bio-)logic chemical gates, as discussed in [35].Here, rather than exploring further that route, we keep the focus on the outlined fundamental bricks and analyze their behavior away from the thermodynamic limit.

The mechanical formulation of chemical kinetics
The analogies between the saturable systems described in the previous section have been developed in [34] restricting to the thermodynamic limit.This is a plausible regime for experiments involving extensive solutions of reactants and, indeed, the thermodynamic limit underlies also standard approaches in early modeling (e.g., classical chemical-reaction kinetics) [24,51,54,55].However, thanks to novel experimental breakthroughs, finally a number of recent studies involves only small numbers of molecules thus questioning the validity of any description in terms of such a large scale limit [56].This broad class of novel experiments includes toggle switches [4,46,57], stochastic bistability [8,22,23], reactions with noise-induced cooperativity [9,10,14] and the whole quorum sensing [12,25,58,59], just to cite a few: interestingly, these novel experiments in small systems have highlighted such complex behaviors which can not be recovered from theories based on the thermodynamic limit.Scope of the current section is thus to extend the mapping described in the previous section in order to include small systems as well.This will be obtained in two steps.First, we need to frame the whole statistical mechanical treatment of chemical kinetics into the mathematical scaffold of classical mechanics; the latter, extensively relying on non linear PDEs techniques, allows for an optimal mathematical control of these systems at finite sizes [60][61][62][63][64].Then, through this route, we extend the statistical mechanical treatment of reaction kinetics even to the case of finite systems and solve for the latter.Finally, we check the overall robustness of our theoretical predictions by recovering these novel outlined phenomena.
In the following we first present the general mapping, then we handle the simpler case N → ∞ to check that the known limit is properly recovered, and finally we address the finite-N case in detail.We emphasize that, in both the regimes, the model is exactly solvable as a consequence of the complete integrability of the PDEs derived for the partition function and the (related) system's free energy.
The differential identities for the partition function.Let us consider the statistical mechanical formulation of a system where binding sites are of two kinds and evaluate the exact partition function and the related free energy in the thermodynamic limit (N → ∞) and in the case of finite volumes (N ∞).The present theory can be straightforwardly extended to account for an arbitrary number of different binding sites [40][41][42] and to address chain reactions (whose implementation can be helpful in several biochemical information processing systems [21,53]).
The general model we consider is described by a Hamiltonian of the form where σ and τ are Ising spins associated to the binding sites of type A and B respectively, J AB , J AA , J BB are the pairwise-coupling constants between sites (of type A and type B, both of type A, both of type B, respectively), and the external fields h A and h B correspond to the chemical potentials associated to the party A and to the party B, respectively.The overall number of sites is N = N A + N B , in such a way that, setting γ = N A /N , we have The relative magnetizations are defined as The partition function for the system can then be written as where we have defined By direct differentiation one can immediately verify that Z satisfies the following set of compatible PDEs Based on formula (27), it can be proven that the solution Z(x A , x B , t A , t AB , t B ) must satisfy the initial condition The free energy F of the system, defined as F = 1 N log Z, consequently fulfils the set of equations with initial condition , where The thermodynamic limit.In this section we evaluate the free energy and the equations of state for the bi-partite system (25) in the thermodynamic limit N → ∞, where the limit is taken keeping the ratio γ = N A /N constant (that is, none of the two parties is negligible w.r.t. the other).Hence, neglecting O(1/N ) terms in equations ( 30), we have the following system of PDEs The equations above are expected to provide an accurate description of the thermodynamic solution within regions of thermodynamic variables {x A , x B , t A , t AB , t B } such that the second order derivatives in equations ( 30) are bounded.We observe that the system of equations ( 32) is a completely integrable set of Hamilton-Jacobi type equations and the general solution can be calculated via the method of characteristics (see e.g.[65]).Hence we find that the general solution can be written as follows where m A and m B are functions of the thermodynamics variables x A , x B , t A , t AB , t B defined by which can also be written as The functions m A,B (x A , x B , t A , t AB , t B ) are obtained extremizing the free energy F (x A , x B , t A , t AB , t B ), i.e.

∂F ∂ m
or, equivalently, where the function w( m A , m B ) is uniquely fixed via the initial condition on m A and m By evaluating equations (34) for the function w given in (35), we obtain the self-consistency equations in the form We observe that the equations of state ( 34) can be interpreted as the hodograph solution (see e.g.[66]) of the following set of (1 + 1)-dimensional equations of hydrodynamic type for the partial magnetizations m A and m The set of equations above is completely integrable as directly follows from the equations (32) with By differentiating equations (36), where magnetizations m A and m B explicitly depend on x and t variables, one can show that all first and second derivatives develop a gradient catastrophe on the hyper-surface The solutions t B , t A , m to the algebraic equation ( 39) are singular points on the space of coupling constants (x A , x B , t A , t B , t AB ) where partial magnetizations develop a classical shock.Such shocks, as they are known in hyperbolic wave theory, are associated to phase transitions in statistical mechanics [41,60,[62][63][64].The equation (39) and the self-consistency equations (36) imply that the point of vanishing magnetization m A = m B = 0 develops a gradient catastrophe if When addressing the comparison of these outcomes with recent experimental results later, we will be focusing on models with no intra-party interactions (zero coupling between spins of the same type), i.e. t A = t B ≡ 0: this means that binding sites of the same type neither cooperate nor compete, while interactions between binding sites of different nature will be retained and can be both cooperative or competitive.In this case the critical point will simplify into By plugging the expression for w appearing in Eq. ( 35) into the equation ( 33), the required solution F reads as follows Recalling (11), the relation between the average magnetizations m A and m B and the variables Y A , Y B is given by Hereafter, in order to lighten the notation, the bracket for Y A,B shall be dropped.Before proceeding it is worth recalling that the free energy can be decomposed into F (β) = −βU + S, where the energetic and the entropic contribution are highlighted.In particular, the former takes the form Following the minimum energy principle we look for the states that minimize U.In particular, as long as J ≥ 0 and h A h B > 0, these states are given by Otherwise stated, if both the types of binding sites display a positive affinity with the ligand and reciprocal cooperativity is either absent or positive, then, according to purely energetic prescriptions, both the types of binding sites tend (not) to bind the ligand if its concentration α is (smaller) larger than α 0 .
Let us now consider the entropic contribution.
To see the effect of the maximum entropy principle at work instead, we differentiate S with respect to both Y A and Y B , and impose the maximum condition, getting Y A = Y B = 1/2: as expected, the states favored by the entropic term are those corresponding to half saturation for both binding sites (the most disordered states available to the system).
In general, when both entropic and energetic contributions are considered, the equilibrium states stem from an interplay between the two (as the various parameters are tuned, i.e. the noise level β, the relative size γ of the two parties, the half-saturation reference α 0 , and the strength of the reciprocal coupling J A , J B , J AB ).
As anticipated, the maximization of F allows accounting for both extremization and this leads to (see the saturation curves in Eq.s (( 36)) and the mapping in Eq. ( 43)) These theoretical behaviors are successfully used to fit experimental data from positive-cooperative and negativecooperative systems (see Fig. 4) and they are also tentatively tested on ultra-sensitive kinetics (see Fig. 5), although ultra-sensitivity requires particular care because of the breakdown of the self-averaging property of the saturation function.
In fact, within the statistical mechanical framework, we know that -away from phase transitions-the magnetization is a self-averaging order parameter, that is, when N → ∞ the distribution of the magnetization becomes delta-peaked on its thermodynamic average value, namely, lim N →∞ P (m) → δ(m− m ), while when N is finite this distribution is only a Gaussian (still centered on m ) whose variance vanishes as N grows scaling as 1/ √ N .However, when J → J c = 1, namely when the system exhibits a phase-transition (or when a shock develops in the language of non-linear waves) the variance diverges as N → ∞.In practice, this means that, away from phase transitions, lim N →∞ ( m 2 − m 2 ) = 0 but, exactly on the critical point (i.e.J = J c and h = 0), this relation does not hold any longer and fluctuations in the order parameter grows indefinitely with the size.Analogous arguments apply to Y as well, since Y and m are linearly related.Then, the criticality (whose signatures emerge even in real systems as N gets large, α = α 0 and J → J c ) has important consequences also from a practical perspective.For any J > J c , as larger and larger systems are considered, the saturation function Y (α) displays a discontinuity at α = α 0 : at that point the system jumps from a (almost) empty state to a (almost) fully-occupied state (the jump is meant as a real discontinuity of Y (α) at α = α 0 only for N → ∞ and the degree of occupation of the two extremal states depends on the level of noise).The lack of smoothness in the function Y (α) makes the application of regression techniques awkward.Thus, despite from a biochemical perspective ultra-sensitive reactions are only particularly strong cooperative reactions, from a mathematical perspective these two types of reactions are actually very different.In the next subsection we will develop the small N theory that offers a more refined level of description for these critical systems, and, in particular, we will succeed in evaluating the "jump" that Y (α) experiences at α = α 0 at finite volume N using shock theory: this will result in a practical instrument that can be used in modern experiments on small systems reaction kinetics.36).The best fit coefficients are J best-fit ≈ 0.45 ± 0.06, 0.24 ± 0.04, 0.27 ± 0.04, giving nH ≈ 1.82 ± 0.20, 1.32 ± 0.07, 1.37 ± 0.08 (see Eq. ( 19)) to be compared with the values n lit H ≈ 1.83 ± 0.28, 1.41 ± 0.10, 1.34 ± 0.16 obtained in [67] through a standard Hill fit.In the right panel we consider three sets of data taken from [68,69].Data ( ) taken from the work by Garnier et al. (see Fig. 4 in [68]) represent the fractional saturation of NAD-dependent glutamate dehydrogenase (GDH) from Laccaria bicolor with glutamate as substrate at pH 7.4.These experimental data are then fitted according to Eq. ( 36).The best fit coefficient is J best-fit ≈ −1.6 ± 0.5 giving nH ≈ 0.38 ± 0.09 to be compared with the value n lit H ≈ 0.3.Data (• and ) taken from the work by Glover et al. (see Fig. 1a and 1e in [69]) represent the amino-acid uptake in B. subtilis strain NP1 against the reciprocal of the concentration of the transport substrate given by L-Arginine and L-Phenylanine, respectively.These experimental data are then fitted according to Eq. ( 36).The best fit coefficients are J best-fit ≈ −1.55 ± 0.40, −2.1 ± 0.25 giving nH ≈ 0.39 ± 0.06, 0.32 ± 0.03 to be compared with the value n lit H ≈ 0.37, 0.31.Notice that, in general, experimental data are reported as a function of α (which is the convention adopted here as well) and different systems display a different value of α0.
Finite-size solution.Here we assume that N is finite and fixed and we study the exact solution of the system (28) with the initial condition (29).Let us remark that the initial condition (29) can be equivalently written as In particular, the finite-size magnetizations take the following explicit expressions Notice that, by fixing t A = t, t AB = t B = 0, x B = 0 and N A = N , we recover the finite-size solution to the Curie-Weiss model, as expected.Also, it is immediate to recognize that exp[Λ k,l (x A , x B , t A , t B , t AB )]/Z plays the role of the finite-size canonical probability for the configuration with k and l bounded sites out of N A and of N B , respectively (corresponding to a magnetizations 2k − N A and 2l − N B ). Fitting Eq.s ( 50) on ultra-sensitive kinetics (as shown in Fig. 5) produces sensibly better results than using the infinite limit counterpart coded by Eq.s 47. ].We report three distinct fits: the dark curve refers to Eq. ( 17), strictly valid in the thermodynamic limit, with relatively small coupling constant (i.e., J best-fit < 1), hence recovering a cooperative system in the thermodynamic limit; the steep curve refers again to Eq. ( 17) where the coupling constant is large (i.e., J best-fit > 1) in such a way that the response of the system gets discontinuous; the bright line refers to Eq. ( 50) (where only one party is retained) valid for finite-size systems and the coupling constant is large (i.e., J best-fit > 1), hence recovering an ultra-sensitive system of finite size.Notice that in the latter case, given the finiteness of the system, the discontinuity is smoothened.Indeed, in the language of statistical mechanics, a genuine first order phase transition, is truly captured only in the thermodynamic limit: this mirrors an infinite cooperativity in the biochemical counterpart as it returns nH → ∞.The relative goodness of the fits are R 2 ≈ 0.85, R 2 ≈ 0.94, and R 2 ≈ 0.99, respectively, confirming an ultra-sensitive behavior.Right: Finite size scaling of the saturation function at different sizes (i.e N = 2, 10, 100 and N → ∞ values are shown) gives useful information on the goodness of the numerical implementation of our approach as, for small systems of concrete interest (e.g.N < 50), the presented procedure is already feasible on standard machines, while for larger N (e.g.N > O( 102 )) the finite-size solution already collapses to its thermodynamic asymptotics.

Implications of the theory for small systems
The successful comparison between the experimental saturation plots and the predictions from self-consistencies obtained in the previous section provides a sound check of the analogy developed.Now, we aim to go further to recover novel emerging phenomena, whose experimental evidence is well documented.In particular, we now focus on cooperative-like and bistable behaviors induced by noise and stochastic effects [4,8,22,23,46,57].Next, this phenomenology is further investigated in its relationship with quorum sensing [12] and we will finally offer a cybernetic interpretation of this kind of process: this will allow us to provide a natural and robust explanation about how, during a biochemical reaction, the signal (carried by the input) is amplified, while an eventual noise becomes spontaneously attenuated.
Noise-induced cooperativity and bistability.It is well-known that cooperative binding in real systems (i.e. at finite sizes) induces a bistable behavior, which we call cooperativity-induced bistability.On the other hand, growing attention has been recently captured by the evidence that bistability can also occur without cooperativity (i.e. it may happen in systems whose binding sites do not interact at all, simply as a consequence of the intrinsic stochasticity, see e.g.[46] and references therein).We call this phenomenon noise-induced bistability.Indeed, these two kinds of bistability display a significant empiric resemblance that makes one speculate that stochasticity can play a a role in the effective cooperativity.However, the two phenomena are conceptually different, as it can be inferred from the present statistical mechanical treatment.The simplest way to do this is by starting from the mono-partite system discussed in Sec. 2 and compare the self-consistent expression (17) with the saturation function of cooperative systems (9).Both expressions are recalled hereafter for clarity: Notice that in Eq. ( 51) the explicit β is retained (i.e., the rescaling βJ → J and βh → h is no longer applied) since here we are focusing on the role of β.Now, it is easy to see that stochastic effects may drift the system away from Michaelis-Menten behavior (toward a cooperative-like one), even if there is no cooperation among binding sites.In fact, by assuming that the coupling is strictly zero, i.e. by setting J ≡ 0 in the self-consistency ( 51), we get whence, recalling that tanh r = 1 − 2/(e 2r + 1), we have The last expression coincides with the Michaelis-Menten law (6) as long as β = 1 and α 0 coincides with the Michaelis constant (defined as the ratio between dissociation and association constants related to the reaction considered, which corresponds to the concentration of the ligand at which the reaction rate is half its maximum value).More generally, letting β vary (hence rising or lowering stochastic effects in the system), we can reshuffle the previous expression as which recovers the Hill law (9) as long as we take K −1 = α β 0 .Focusing on the dependence of Y on the concentration α, we find that β plays the same role as the Hill coefficient: the smaller the stochasticity (i.e. the larger β) the more "cooperative" the system.The reason for this behavior is apparent in the statistical mechanical picture, as low levels of noise make the spins of the system align faster (driven by the external field, rather than by themselves) by reducing overall the fluctuations.Of course, the displayed behavior depends globally, rather than locally, on the system energy, i.e. it may not be ascribed to the mutual interaction of sites, as it is the case for truly cooperative systems.
In Fig. 6 (left panel) we show the qualitative behavior of the saturation function Y as a function of the concentration α and with noise parameter β, according to Eq. ( 55).
Let us now proceed by showing that noise can even induce bistability.To this aim, we need to extend this approach to systems built of by two-parties (namely to the general model coded in Eq. ( 25)); despite being mathematically more involved, the reasoning for multipartite systems goes as much as the same as above.Indeed, profiting the mappings (43), the partial saturation curves may be shown to represent Hill law for each party also in the case of two (or several) interacting sites.
To see this, let us set where α j denotes the concentration of free molecules for the party j = A, B and α j 0 denotes the reference value representing equal probability of being occupied and unoccupied for the party j.  55).At high level of noise the Michaelis-Menten envelop is preserved (the curve at β = 1/2), while for higher values of β -hence for smaller noises -a sigmoidal shape slowly emerges (i.e. for β = 2, 4).(Right Panel) Noise induced bistability.The saturation function of a two-type site reaction YA(α), YB(α) is shown for various values of the noise β, according to Eq. ( 56).Again, while for small values (i.e.β = 1/2) the system behaves as the superposition of two independent Michaelis-Menten reactions, for higher values of β (i.e.β = 4) bistability clearly appears in the system.Here we used αA = αB and α A 0 = α B 0 = 1/2.
Assuming again J = 0 and N A = N B = N/2 for the sake of simplicity, the event of each site being bound (respectively unbound) is independent of the state of the other sites, thus, denoting with σ j i the i th spin (i = 1, . . ., N j ) of the party j, the total concentration factorizes due to the independence of the probabilities (a major advantage of the mean-field approach) as , whence α tot = α A • α B (and the same conclusion may be drawn for the total reference value α tot 0 ).As a consequence of the superposition principle for the external fields and of (13), we then have whence, having set N A = N B , we can put by symmetry (see Eq. ( 13)) Exploiting again the assumption J = 0, that is, in the absence of cooperativity Eq. ( 47) reads as whence, replacing the hyperbolic tangent with its exponential expression, we recover the partial Hill laws for the single parties.Notice that, in order for the analogy with eq. ( 9) to be preserved, the partial Hill coefficients coincides now with 2β, rather than β, for each single party consists of precisely one half of the total number of spins.This is indeed perfectly consistent with our previous remark on the global dependence on the noise β in the case of noise-induced "cooperativity" and is indeed in sharp contrast with the property of the truly cooperative behavior to be invariant with respect to the system size.
In order to graphically illustrate the properties of the partial saturation curves Y j with respect to the noise β, fix e. g. h A = h B = h tot /2, so that α B = α A α B 0 /α A 0 and α tot = (α A ) 2 α B 0 /α A 0 .By fixing a value of α tot (or, which is the same, the values of α j 0 ), we can interpret the α j 's to be relative concentrations with respect to α tot , hence rewriting Y B (α B ) = Y B (α tot /α A ) ≡ Y B (α A ) for α A ∈ (0, 1).This stochastic-driven behavior for two-party systems is shown in Fig. 6 (right panel).The physical reason for this switches lies in the intrinsic ergodicity of any dynamics involving systems with finite size: considering the simplest bistable system, its two free energy minima are separated by a maximum (i.e., a barrier) whose height is N δF and the characteristic time τ for the system to move from one minimum to another typically scales as τ ∝ exp(δF ) thus, solely in the thermodynamic limit, this barrier becomes infinite and the system gets trapped forever into one of these minima (and ergodicity becomes broken): this is clearly discussed from a biochemical perspective in e.g., [16,46] (and references therein).We checked that our theory correctly reproduces these spontaneous switches in time, as reported in Fig. s 7 and 8 (left panel) where we recover the original plots shown [46] using the same parameters.
In our vocabulary, this happens because the asymptotic evaluation of the magnetization already for the simplest bistable mono-partite spin model (i.e., a Curie-Weiss model) leads to the expansion of the form for large N , where ξ is the solution to the self-consistency equation such that the free energy attains its minimum.Away from the thermodynamic limit, keeping only the first order correction to finite volume, we can approximate and η can be interpreted as a stochastic contribution.Indeed, the magnetization m has zero-variance in the limit N → ∞ only (or in the pathological case of zero noise η ≡ 0).Solving the equation ( 59) for ξ and substituting in (58), In this way, four distinct regimes can be outlined: at the beginning both parties display half saturation and negative cooperation, being subject to small, positive fields and, as a result, the output associated to the party experiencing the larger field grows while the other decreases; then, both fields are set to zero and the outputs remains close to 1 and to 0, respectively; next, the coupling between the two parties is reduced and the related output collapse to 1/2; now, even if the coupling is restored to large, negative values the two outputs remain null.This phenomenology recovers the original picture by Gardner et al. (see [6], Fig. 5).Right: First order corrections to the infinite volume expression for the transfer function of the Curie-Weiss model: η encodes for random thermal fluctuations that, in this context, are not washed away and actually break locally (in time) the symmetry between the two free energy minima allowing the system to oscillate between them (the stronger the thermal noise the faster the hopping rate between the two minima).
we obtain (up to O(1/N )) Thus, as shown in Fig. 8 (right panel) a non-vanishing η breaks the symmetry of the isothermal curve for the average magnetization m .In absence of external field, i.e. x = 0, a positive (negative) value of η selects a negative (positive) solution m , hence random fluctuations of η are responsible for possible switches of the magnetization thus effective stochastic bistabilities typical of toggle switches.
To further check that our theory correctly reproduces such a flip-flop (i.e.switch) like behavior, we can study the null-clines of the simplest Langevin dynamics coupled to our system (25), namely where τ is the typical time constant (see Fig. 7), and the randomness tuned by η is standard white noise, namely η A (t)η B (t ) ∝ δ(A − B)δ(t − t ): keeping C A , and C B to label the two reactant's concentrations as in the original paper [22], we can compare the null-clines of our system to those pertaining to the following archetypal toggle where we have kept C A and C B to label the un-normalized reactant's concentrations.The corresponding by now well-known results are illustrated in Fig. 9 (lower panel) for the case n and F B analogously), the system's (in-)stability at equilibria may be checked with the sign of the differential's eigenvalues an equilibrium point being a minimum (i.e.stable) if both the eigenvalues are positive and a maximum (i.e.unstable) if both of them are negative, or otherwise a saddle point (unstable) if they are opposite in sign.62), roughly displaying the same behavior as (60) (to be compared to [22], Fig. 2).Notice that in the present case, the concentrations CA, CB are not normalized as they previously were, so that, after a suitable renormalization, the equilibria belong in fact to region [0, 1] × [0, 1].
Condorcet theory: stochastic signal amplification and noise suppression.The Condorcet theorem was originally stated within a political context and concerns the relative probability of a given group of individuals arriving at a correct (binary, i.e.YES-NO) decision.The assumptions underlying the simplest version of the theorem is that a group wishes to reach a decision by majority vote.One of the two outcomes of the vote is correct, and each voter has an independent probability p of voting for the correct decision.The theorem states that if p > 1/2 (i.e., each voter is more likely to vote correctly), then adding more voters to the pool increases the probability that the majority decision is correct.In the thermodynamic limit, the probability (that the majority votes correctly) approaches 1 as the number of voters diverges.On the other hand, if p < 1/2 (i.e., each voter is more likely not to vote correctly), then adding more voters makes things worse.
Condorcet theorem is therefore a majority rule that can be used as a tool to infer collective reliable predictions about the better of two options and, in recent times, deep analogies between Condorcet theory in social systems and quorum sensing in biological information processing systems have been developed: for instance, many species of bacteria use quorum sensing to coordinate gene expression according to the density of their local population [58], and decisions within the adaptive response of the immune system in mammals requires the implementation of quorum sensing by lymphocytes [12,59] (that globally vote if attacking or resting when a novel antigen is presented to them and dialogue among themselves respecting the rules of biochemical reactions).In any case, this distributed decisional capability should emerge as a natural consequence of the underlying reaction kinetics of all the involved agents, thus we should be able to understand its genesis from our general treatment.
Note at first that this mechanism is also consistent with the biologically-adapted [18,34] electronic picture of signal amplification and noise suppression (where, in this analogy, the meaning of the signal is the correct vote and that of noise is the wrong one) as, actually, this is not too far from the Shannon coding theorem in Cybernetics: roughly speaking the latter states that, even if there is a huge noise on a cable (i.e. in the system), it is enough that the probability to send successfully a message through it remains strictly greater than one-half to be sure that the the whole information that crossed the system cannot get lost (with large enough trial samples N → ∞).Indeed in all these saturable systems (as also voter-like models a' la Condorcet are saturable systems by definition if the vote is binary) when collective features are at work, the correct output "yes" emerges neatly even at a relatively mild input, and the wrong output "no" is cleverly avoided even in the presence of significant noise, as shown in Fig. 10.
This kind of phenomenology is naturally captured within our chemical kinetics framework: in fact, in our system the stimulus (i.e., the logarithm of the ligand concentration) is provided by the field h and, according to its effective magnitude, it is expected to return a bound ("yes") or a not-bound ("no") state for the ligands such that, if the stimulus is poor (i.e., chemical noise), it is suppressed in the output, while if the stimulus is relatively consistent (i.e., chemical signal), it gets amplified in the output, as discussed in detail in Fig. 10.Clearly the larger the coupling value J > 0 (or, alternatively, the higher the Hill coefficient of the reaction), the stronger the resulting amplification (see also Eq.s (23)).This point is further deepened hereafter.If we look at the expression for the system's energy (Eq.( 44)), a key point is that when J > 0, the contribution 2JY (1 − Y ) provides a boost (i.e. a gain) for the system's response, further stabilizing the states Y = 0 (for low levels of input -i.e.h < 0-, that is interpreted as chemical noise and it is thus suppressed) and Y = 1 (for high level of input -i.e.h > 0-, that is interpreted as a chemical signal and it is thus amplified).The usefulness of this Condorcet-like mechanism at work in biochemistry becomes evident already in the paradigmatic case of hemoglobin (a cooperative protein responsible for oxygen transport in tissues): when in the lungs (rich of oxygen), hemoglobin uses cooperativity to bind to as much oxygen as possible; when in the tissues (poor of oxygen) it uses cooperativity to get rid of it, thus releasing oxygen in the tissue.

Discussion
In this manuscript we deepened our translation of biochemical kinetics into a statistical mechanical scaffold started in [34,35] and that allows a straightforward cybernetic interpretation of these phenomenologies.The final goal is to contribute in the quantitative understanding of the emergent computational properties possibly shown by large networks of biochemical reactions regarding cell signalling.The route we paved is based on the one-to-one, robust and sharp, structural and behavioral analogy between the response functions in biochemical reactions (i.e.saturation curves), the response functions of mean-field models in statistical mechanics (i.e.self-consistencies), and the response functions of key amplifiers in electronics (i.e.transfer functions).In particular, the response functions of cooperative (anticooperative) systems match those pertaining to mean-field ferromagnets (antiferromagnets), which, in turn, overlap those characterizing an operational amplifier (flip-flop).We decided to keep the discussion at the mean-field level and in the canonical ensemble because it is within these limits that other theories regarding information processing systems (e.g., neural networks in Artificial Intelligence [28]) have been developed in the past, and it is by comparing our results to these theories that we aim to learn how information is treated during biochemical reactions.Notice that, from a statistical mechanical perspective, neural networks are particular types of spin-glasses, namely tricky realizations of ferromagnets and antiferromagnets broadly interacting, while from an electronic perspective, neural networks are realized by suitably combining large numbers of amplifiers and flip-flops, thus the first steps in this structural equivalence should focus on these basic ingredients, that have been indeed the subject of the present paper.Moreover, many of the biochemical reactions of current interest involve small numbers of agents and this makes theories in the thermodynamic limit [56] too coarse for the purpose of tackling the proposed analogy in full generality.Therefore, in this manuscript we tried to fix this point by developing a proper mathematical technique able to account even for small-sized systems.Summarizing the whole procedure, our route first requires mapping the original biochemical problem into a Hamiltonian formulation; the resulting model can then be embedded into a statistical mechanical framework.Next, the solution for this kind of systems is attainable by noting that their related free energies obey Hamilton-Jacobi type equations in the space of their parameters (ultimately mapping a problem in biochemistry into a problem of analytical mechanics): relying on the complete integrability of the system of multidimensional PDEs, we solve the general scenario in both the thermodynamic limit and the finite size case [65].In particular, we observed that the multidimensional equations of state can be constructed via the hodograph equations that in turn provide the solution to a system of hydrodynamic type [66].Following this route we are able to provide a theoretical description of several complex phenomena stemming from finite-size effects.This phenomenology is rather broad and includes the effective cooperativity induced by stochasticity, as well as an enhanced bistability when dealing with toggles.Crucially, our procedure naturally allows a further interpretation of the response functions of these biochemical systems in terms of cybernetics, that constitutes a novel and transparent way to analyze how information is handled during these reactions: once shown the one-to-one behavioral correspondence between saturation functions in chemical kinetics, self-consistencies in statistical mechanics and transfer function in electronics (these are all identical saturable response functions from a cybernetic perspective), we related Condorcet phenomenology with signal amplification and noise suppression and framed it into the elementary scenario of ferromagnetic gain.We tested our theory both against classical experiments (i.e., in the large N limit), recovering all the main equations of reaction kinetics as suitable limits (i.e.Micaelis-Menten, Hill, Adair equations) as well as against novel experiments involving small system sizes, recovering the expected phenomenology [6,16,46].Further developments of the theory now should proceed in three ways: on one side, still keeping the Maxwell-Boltzmann prescription, we can combine small (bio-chemical) circuits together in order to analyze information processing in larger chemical networks.On the other side, efforts are still needed to enlarge this scheme in order to apply to general out-of-equilibrium regimes (for instance with time-dependent field variations).Finally the above procedure can be further enriched by working in synergy with other well developed (and possibly alternative) mathematical methods, especially those already thermodynamically oriented -i.e., thougth to deal with the complexity of evaluating the partition function at finite size-as, for instance, those geometrically-oriented reported in the review by Ruppeiner [70].

Figure 1 :
Figure 1: Behavioral and formal analogies highlighted and exploited in this work.Left panel: the whole logical procedure to understand information processing during chemical reactions is schematically shown.(a) We deal with a problem in biochemistry.(b) Then we model it into a statistical mechanical framework, that (c) it is then solved via techniques of analytical mechanics and (d) whose results are finally interpreted in cybernetical/electronical terms.These findings can be further translated back into the biochemical scenario, whose modus operandi has now become transparent.Right panel: Relevant behavioral analogies in the response of these saturable different systems lie at the basis of the structural equivalence we use and celebrated examples are reported.(a) The sigmoidal shape of the transfer function of an operational amplifier, where the response is the output voltage while the stimulus is the input voltage.(b) The sigmoidal shape of the self-consistency of a ferromagnet, where the response is the magnetization and the stimulus is the external magnetic field.(c) The sigmoidal response of the activation function of a neuron where the output is the action potential voltage while input is conveyed by all the afferent electric synaptic currents.(d) The sigmoidal shape of the saturation curve of a cooperative chemical reaction, where the output is the fraction of bound sites, while the input is the (log-)concentration of the ligand.

Figure 2 :
Figure2: This figure summarizes the structural analogy between biochemical and electronic information processing.In the top-line three different proteins are shown: from left to right, the first (mitogen-activated protein kinase 14) obeys cooperative kinetics, the second (calmodulin dependent protein kinase 2) ultra-sensitive kinetics, while the last (synaptic glutamate receptor) fulfills anti-cooperative kinetics.In the second row the saturation curves for these proteins are shown: α is the concentration of the free ligand needed for their binding and Y is the fraction of occupied binding sites.Symbols with the relative error-bars stand for real data taken from[20,38,39] and lines are best fits performed through the analytical expression in Eq. (17), where the best-fit coefficient J corresponds to an ergodic ferromagnet (left), a low-temperature ferromagnet (center) and an anti-ferromagnet (right).In the third row three circuits are shown.From left to right, an operational-amplifier, an analogto-digital converter and a flip-flop, while in the fourth row their transfer functions are presented to highlight the behavioral analogy with the proteins.The latter shows the response of the system (output voltage) as a function of the input of the system (input voltage).Following columns instead of rows, the first column is due to cooperative behavior, the second to ultrasensitive behavior and the last to anticooperative behavior.The authors are grateful to Nature Publishing Group for the permission to reproduce this figure from[34].
by means of illustrative examples.Once wrote the general theory, at first we successfully compare the expected behavior of our theoretical saturation function with experimental results for several examples of cooperative, anticooperative and ultrasensitive kinetics of

Figure 3 :
Figure 3: Schematic representation of the processes considered.The cooperative case is shown in the left side, while the anti-cooperative case is shown in the right side.Two different binding sites, corresponding to two spins belonging to different parties, are shown in different colors.The cooperative case is characterized by a Hill coefficient nH larger than 1, corresponding to a positive coupling J between spins.The anti-cooperative case is characterized by a Hill coefficient nH smaller than 1, corresponding to a negative coupling J between spins.The input of the system is provided by the concentration of free ligand α, corresponding to the magnetic field h (more precisely, h = 1/2 log(α/α0), where α0 is the half-saturation concentration, in such a way that a negative field is equivalent to a concentration smaller than α0 and vice versa).The output is provided by the saturation Y , corresponding to the magnetization m (more precisely, m = 2Y − 1, in such a way that the mean performed over the magnetizations of the two parties is equivalent to the sum of the saturations pertaining to the two binding sites minus 1).The plots in the bottom represent the evolution of the output as the input is varied.In particular, the curves shown here are obtained as the solution of Eq. (17) for J = 0.5 and J = −1.2,respectively.

Figure 4 :
Figure 4: Examples of accordance between our theory and experimental results for positive cooperativity (left panel) and negative cooperativity (right panel).Left panel: data taken from the work by Watson et al. (see Fig. 5 in [67]) representing the fractional saturation of immobilized glucocorticoid receptor (GR) surfaces at 35 • C as the concentration of GR DNA-binding domain (DBD) is varied; three series are considered: WT-Gilz (♦), A477T-Pal-R ( ), and A477T-Gilz ( ).These experimental data are then fitted according to Eq. (36).The best fit coefficients are J best-fit ≈ 0.45 ± 0.06, 0.24 ± 0.04, 0.27 ± 0.04, giving nH ≈ 1.82 ± 0.20, 1.32 ± 0.07, 1.37 ± 0.08 (see Eq. (19)) to be compared with the values n lit H ≈ 1.83 ± 0.28, 1.41 ± 0.10, 1.34 ± 0.16 obtained in[67] through a standard Hill fit.In the right panel we consider three sets of data taken from[68,69].Data ( ) taken from the work by Garnier et al. (see Fig.4in[68]) represent the fractional saturation of NAD-dependent glutamate dehydrogenase (GDH) from Laccaria bicolor with glutamate as substrate at pH 7.4.These experimental data are then fitted according to Eq. (36).The best fit coefficient is J best-fit ≈ −1.6 ± 0.5 giving nH ≈ 0.38 ± 0.09 to be compared with the value n lit H ≈ 0.3.Data (• and ) taken from the work by Glover et al. (see Fig.1aand 1e in[69]) represent the amino-acid uptake in B. subtilis strain NP1 against the reciprocal of the concentration of the transport substrate given by L-Arginine and L-Phenylanine, respectively.These experimental data are then fitted according to Eq. (36).The best fit coefficients are J best-fit ≈ −1.55 ± 0.40, −2.1 ± 0.25 giving nH ≈ 0.39 ± 0.06, 0.32 ± 0.03 to be compared with the value n lit H ≈ 0.37, 0.31.Notice that, in general, experimental data are reported as a function of α (which is the convention adopted here as well) and different systems display a different value of α0.

Figure 5 :
Figure 5: Left: Example of accordance between our theory and experimental results for ultra-sensitive systems.Data ( ) are taken from the work by Bradshaw et al. (see Fig. 3 in [20]) representing CaMKII auto-phosphorylation level at equilibrium versus [Ca 2+].We report three distinct fits: the dark curve refers to Eq. (17), strictly valid in the thermodynamic limit, with relatively small coupling constant (i.e., J best-fit < 1), hence recovering a cooperative system in the thermodynamic limit; the steep curve refers again to Eq. (17) where the coupling constant is large (i.e., J best-fit > 1) in such a way that the response of the system gets discontinuous; the bright line refers to Eq. (50) (where only one party is retained) valid for finite-size systems and the coupling constant is large (i.e., J best-fit > 1), hence recovering an ultra-sensitive system of finite size.Notice that in the latter case, given the finiteness of the system, the discontinuity is smoothened.Indeed, in the language of statistical mechanics, a genuine first order phase transition, is truly captured only in the thermodynamic limit: this mirrors an infinite cooperativity in the biochemical counterpart as it returns nH → ∞.The relative goodness of the fits are R 2 ≈ 0.85, R 2 ≈ 0.94, and R 2 ≈ 0.99, respectively, confirming an ultra-sensitive behavior.Right: Finite size scaling of the saturation function at different sizes (i.e N = 2, 10, 100 and N → ∞ values are shown) gives useful information on the goodness of the numerical implementation of our approach as, for small systems of concrete interest (e.g.N < 50), the presented procedure is already feasible on standard machines, while for larger N (e.g.N > O(10 2 )) the finite-size solution already collapses to its thermodynamic asymptotics.

Figure 6 :
Figure 6: (Left Panel) Noise induced cooperativity.The saturation function of a single-type site reaction Y (α) is shown for various values of the noise β, according to Eq. (55).At high level of noise the Michaelis-Menten envelop is preserved (the curve at β = 1/2), while for higher values of β -hence for smaller noises -a sigmoidal shape slowly emerges (i.e. for β = 2, 4).(Right Panel) Noise induced bistability.The saturation function of a two-type site reaction YA(α), YB(α) is shown for various values of the noise β, according to Eq. (56).Again, while for small values (i.e.β = 1/2) the system behaves as the superposition of two independent Michaelis-Menten reactions, for higher values of β (i.e.β = 4) bistability clearly appears in the system.Here we used αA = αB and α A 0 = α B 0 = 1/2.

Figure 7 :
Figure 7: Upper panel: example of historical series for a toggle switch we simulated through a bipartite spin system with overall N = 40, NA = NB = 20 and JAB = −2.8,JA = JB = 0.The behavior of this stochastic flip-flop is as much the same as that of classical toggle switches, as it is evident comparing this plot with Fig. 3 of [46] or Fig 3(a) of [16].Lower panels: the switching time for various couplings' strength is analyzed.The complementary of the cumulative distribution Cum(τ ) is shown in the left for JAB = −3.0,−2.8, −2.7, −2.5 (ordered from the least to the most steep curve).The scaling law τ = aJ exp(bJ N ) is checked in the right panel, where, again different coupling strengths are considered [JAB = −3.0(), −2.8( ), −2.7( ), −2.5(•)]; to be compared to [16] Fig. 3b.

Figure 8 :
Figure 8: Left: Temporal behavior of the outputs YA and YB for the two parties.Along the time, the parameters J, hA and hB are varied as specified in the figure.In this way, four distinct regimes can be outlined: at the beginning both parties display half saturation and negative cooperation, being subject to small, positive fields and, as a result, the output associated to the party experiencing the larger field grows while the other decreases; then, both fields are set to zero and the outputs remains close to 1 and to 0, respectively; next, the coupling between the two parties is reduced and the related output collapse to 1/2; now, even if the coupling is restored to large, negative values the two outputs remain null.This phenomenology recovers the original picture by Gardner et al. (see[6], Fig.5).Right: First order corrections to the infinite volume expression for the transfer function of the Curie-Weiss model: η encodes for random thermal fluctuations that, in this context, are not washed away and actually break locally (in time) the symmetry between the two free energy minima allowing the system to oscillate between them (the stronger the thermal noise the faster the hopping rate between the two minima).

Figure 9 :
Figure 9: Upper panel: Null-clines of Langevin equations system(60), displaying equilibria and their behavior, with varying J. Lower panel: Null-clines of the system (62), roughly displaying the same behavior as (60) (to be compared to[22], Fig.2).Notice that in the present case, the concentrations CA, CB are not normalized as they previously were, so that, after a suitable renormalization, the equilibria belong in fact to region [0, 1] × [0, 1].

Figure 10 :
Figure 10: Left panel: Transfer function for cooperative systems (i.e. chemical cooperative kinetics, physical ferromagnets, electronic operational amplifiers).Curves pertaining to different interaction strength among the system's units (namely pertaining to a different Hill coefficient, or a different ferromagnetic coupling, or a different gain, according to the context considered) are shown in different colors.The dashed line highlights the linear behavior expected in the region of small signal, namely for a perfect cable.As the interaction strength gets larger we eventually get a switch.The Condorcet scenario results in both signal amplification at high ligand concentration as well as noise suppression at low ligand concentration.Right panel: Amplification scheme.In the main plot the linear regime is highlighted.When the input (i.e., the external field h, left inset) is varied within this region the resulting output (i.e., the magnetization m, right inset) varies according to a linear relation.Outside the linear regime the input-output relation can exhibit distortions.