Complex Reaction Kinetics in Chemistry: A unified picture suggested by Mechanics in Physics

Complex biochemical pathways or regulatory enzyme kinetics can be reduced to chains of elementary reactions, which can be described in terms of chemical kinetics. This discipline provides a set of tools for quantifying and understanding the dialogue between reactants, whose framing into a solid and consistent mathematical description is of pivotal importance in the growing field of biotechnology. Among the elementary reactions so far extensively investigated, we recall the socalled Michaelis-Menten scheme and the Hill positive-cooperative kinetics, which apply to molecular binding and are characterized by the absence and the presence, respectively, of cooperative interactions between binding sites, giving rise to qualitative different phenomenologies. However, there is evidence of reactions displaying a more complex, and by far less understood, pattern: these follow the positive-cooperative scenario at small substrate concentration, yet negative-cooperative effects emerge and get stronger as the substrate concentration is increased. In this paper we analyze the structural analogy between the mathematical backbone of (classical) reaction kinetics in Chemistry and that of (classical) mechanics in Physics: techniques and results from the latter shall be used to infer properties on the former.


A. The Chemical Kinetics background
The mathematical models that describe reaction kinetics provide chemists and chemical engineers with tools to better understand, depict and possibly control a broad range of chemical processes (see e.g., [27,48]). These include applications to pharmacology, environmental pollution monitoring, food industry, etc. In particular, biological systems are often characterized by complex chemical pathways whose modeling is rather challenging and can not be recast in standard schemes [2,3,9,10,20,23,28,36,38,42,43,55,56] (see also [57,59,60,64] for a different perspective). In general, one tries to split such sophisticated systems into a set of elementary constituents, in mutual interaction, and for which a clear formalization is available [4,6,39,44,58,62].
In this context, one of the best consolidated, elementary scheme is given by the Michaelis-Menten law. This was originally introduced by Leonor Michaelis and Maud Menten to describe enzyme kinetics and can be applied to systems made of two reactants, say A (the binding molecule or, more generally, the binding sites of a molecule) and B (the free ligand, i.e., the substrate), which can bind (and unbind) to form the product AB. If we call S the concentration of free ligand, Y the saturation function (or fractional occupancy), namely the fraction of bound molecules (Y ∈ [0, 1]), and, accordingly, 1 − Y the fraction of the unbound molecules, under proper assumptions, one can write where k is the proportionality constant between response and occupancy (otherwise stated, it is the ratio between the dissociation and the association constants). In particular, as standard, it is assumed that the (a) the reaction is in a steady state, with the product being formed and consumed at the same rate, (b) the free ligand concentration is in large excess over that of the binding molecules in such a way that it can be considered as constant along the reaction, (c) all the binding molecules are equivalent and independent. Also, the derivation of the Michaelis-Menten law is based on the law of mass action.
Reshuffling the previous equation we get Y = S/(S + k) which allows stating that k is the concentration of free ligand at which 50% of the binding sites are occupied (that is, when S = k, then Y = 1/2). Thus, denoting with S 0 the half-saturation ligand concentration, we get This equation represents a rectangular hyperbola with horizontal asymptote corresponding to full saturation, that is Y = 1; this is the typical outcome expected for systems where no interaction between binding sites is at work [63]. This model immediately settled down as the paradigm for Chemical Kinetics, somehow similarly to the perfect gas model (where atoms, or molecules -collisions apart -do not interact) of the Kinetic Theory in the early Statistical Physics [49]. Nevertheless, deviations from this behaviour were not late to arrive: the most common phenomenon was the occurrence of a positive cooperation among the binding sites of a multi-site molecule. Actually, many polymers and proteins exhibit cooperativity, meaning that the ligand binds in a non-independent way: if, upon a ligand binding, the probability of further binding (by other ligands) is enhanced, the system is said to show positive cooperativity.
To fix ideas, let us make a practical example and let us consider the case of a well-known protein, i.e., the hemoglobin. This is responsible for oxygen transport throughout the body and it ultimately allows cellular respiration. Such features stem from hemoglobin's ability to bind (and to dislodge as needed) up to four molecules of oxygen in a nonindependent way: if one of the four sites has captured an oxygen molecule, then the probability that the remaining three sites will capture further oxygen increases, and vice versa. As a result, if the protein is in an environment rich of oxygen (e.g., in the lungs), it readily binds up to four molecules of oxygen, and, as much readily, it gets rid of them when crossing an oxygen-deficient environment. To study quantitatively its behaviour one typically measures its characteristic input-output relation. This can be achieved by considering a set of M elementary experiments where these proteins, in the same amount for each experiment, are prepared in a baker and allowed to bind oxygen, which is supplied at different concentrations S i for different experiments (e.g., S 1 < S 2 < ... < S M ). We can then construct a Cartesian plane, where on the abscissas we set the concentration of the ligand S (oxygen in this case, i.e. the input) while on the y-axes we put the fraction of protein bound sites Y (the saturation function, i.e., the output). In this way, for each experiment, once reached the chemical equilibrium, we get a saturation level and we can draw a point in the considered Cartesian plane; interpolating between all the points a sigmoidal curve will emerge (see Fig. 1). Archibald V. Hill formulated a description for the behavior of Y with respect to S: the so-called Hill equation empirically describes the fraction of molecules binding sites, occupied by the ligand, as a function of the ligand concentration [7,8,21,65]. This equation generalizes the Michaelis-Menten law (2) and reads as where n H is referred to as Hill coefficient and can be interpreted as the effective number of binding sites that are interacting each other. This number can be measured as the slope of the curve log[Y /(1−Y )] versus log(S), calculated at the half-saturation point. Of course, if n H = 1 there is no cooperation at all and each binding site acts independently of the others (and, consistently, Michaelis-Menten kinetics is restored), viceversa, if n H > 1, the reaction is said to be cooperative (just like in hemoglobin), and if n H 1 the cooperation among binding sites is so strong that the sigmoid becomes close to a step function and the kinetics is named ultra-sensitive.
The Michaelis-Menten law together with the extension by Hill, provided a good description for a bulk of chemical reactions, however, things were not perfect yet. For instance, some yeast's proteins (e.g., the Glyceraldehyde 3-Phosphate Dehydrogenase [24]) produced novel (mild) deviations from the Hill curve: for these enzymes, the cooperativity of their binding sites decreases while increasing the ligand concentration. The following work by Daniel E. Koshland allowed understanding this kind of phenomenology by further enlarging the theoretical framework through the introduction of the concept of negative cooperativity. In fact, in the previous example, beyond the positive cooperation between the binding sites there are also negative-cooperative effects underlying. Their effective action is to diminish the overall binding capabilities of the enzyme and thus to reduce the magnitude of its Hill coefficient. For the latter we report two fits: Dashed line is the result obtained by constraining the system to be cooperative but not-ultrasensitive (that is, J ≤ 1), while solid line is the best fit (without constraints) which yields to J ∼ 1.1, hence a "first order phase transition" in the language of statistical mechanics. The relative goodness of the fits are R 2 coop ∼ 0.85 and R 2 ultra ∼ 0.94, confirming an ultra-sensitive behavior. The tables in the bottom present the value of J derived from the best fit and the resulting nH ; the estimate of the Hill coefficient taken from the literature is also shown for comparison. This figure was presented in [7].

B. The Mechanics background
The progressive enlargement of a theoretical scaffold to fit the always increasing amount of evidences is a common feature in the historical development of scientific disciplines [34,41]. This is the case also for Mechanics and, as we will see, the analogy with Chemical Kinetics goes far beyond this feature. Beyond Kinematics, which describes the motion of systems without considering their mass or the forces that caused the motion, in the seventeen-th century Newton gave a sharp description of Mechanics, in the form of laws describing how masses dynamically respond when stimulated by an external force (or moment). Here, the input is the force while the output is the motion of the body. The Newtonian dynamics has been ruling for centuries and, in fact, it was so well-consolidated that scientists, among which Giuseppe L. Lagrange, William R. Hamilton and Carl G.J. Jacobi, later reformulated the entire theory in a powerful and elegant variational flavor. The theory was overall brilliant to explain the perceivable reality, but with exceptions emerging in the limit of too little or too fast. We will focus on the latter. In the Newtonian world, if an applied force is kept constant over a mass, this will constantly accelerate, eventually reaching diverging velocities. This was perfectly consistent with the general credo that the speed of light were infinite. However, this postulate broke down in 1887 when the famous experiment by Albert A. Michelson and Edward Morley proved that such a velocity is actually finite. The next years were dense of novel approaches and ideas by many scientists, as Hendrik Lorentz and Hermann Minkowski, and culminated with the special relativity by Albert Einstein in 1905. According to this theory, no mass can exist whose velocity may diverge, the limiting speed being the speed of light. The classical Hamilton-Jacobi equations and Galilean transformations left the place to the Klein-Gordon formulation and Lorentz covariances and contravariances (the natural metric being Minkowskian) [46]. Clearly, classical mechanics was still a good reference framework for the vast majority of the data collected (much like the positive cooperativity accounted for the bulk of the empirical data in the chemical counterpart), however there were rare phenomena (e.g., a muon decay in atmosphere [18]) that required a broader scaffold which, in the opportune limits, could recover the classical one.
Despite this historical connection between Chemical Kinetics and Classical Mechanics may look weird at a first glance, as we will prove, there is a formal analogy between their mathematical representations. In the next section we will summarize the main results concerning the analogy at the classical level. More sharply, the saturation plot of classical (positive cooperative) chemical kinetics (namely the input-output relation between the saturation function and the concentration of the substrate) can be derived by a minimum action principle that is the same that holds in classical mechanics, when describing a mass motion in the Hamilton-Jacobi framework. In this parallelism, the saturation function in Chemistry plays as the velocity in Physics: thus, exactly as happens in special relativity, the velocity of the mass is bounded (by definition, the saturation function can not exceed one). Indeed, we can follow this mathematical equivalence and verify that there is actually a natural broader formulation for chemical kinetics that is exactly through the Klein-Gordon setting (rather than its classical Hamilton-Jacobi counterpart) and the theory as a whole is Lorentz-invariant. Remarkably, when read with chemical glasses, this extended relativistic setting allows for the anti-cooperative corrections that Koshland revealed in the study of the yeast enzymes, resulting in a complex mixture of positive and negative cooperation among binding sites.

II. THE STANDARD MATHEMATICAL SCAFFOLD FOR CLASSICAL COOPERATIVITY
As anticipated in Sec. I A, cooperativity is a widespread phenomenon in Chemistry and its underlying mechanisms can be multiple: for example, if the adjacent binding sites of a protein can accommodate charged ions, the attraction/repulsion between the ions themselves may result in a positive/negative kinetics; in most common cases, the bonds with the substrate modify the protein conformational structure, by influencing possible further links in an allosteric way [6,52]. Whatever the origin, cooperativity in Chemistry is a typical emergent property that directly relates the microscopic description of a system at the single binding-site level, with the macroscopic properties shown by its constituent molecules, cells, and organisms, thus the use of Statistical Physics for its investigation appears quite natural [7,63]. Usually, in Statistical Physics one is provided with an (inverse) temperature β, and with an Hamiltonian (i.e., a cost function) H(σ, J, h) describing the model at the microscopic level, namely in terms of elementary variables σ i , σ j , couplings among elementary variables J ij , and external fields h i acting over these. The goal is to obtain the free energy A(β, J, h) of the model, from which the average value of the macroscopic observables can be directly derived [63].

A. Formulation of the problem: the thermodynamical free energy
In the following we summarize the minimal assumptions needed when modelling chemical kinetics from the Statistical Physics perspective; for a more extensive treatment of this kind of modelling we refer to [5][6][7][8]63], while for a rigorous explanation of the underlying equivalence between Statistical Mechanics and Analytical Mechanics we refer to the seminal works by Francesco Guerra [32], dealing with the Sherrington-Kirkpatrick model (and then deepened in e.g., [14,16,17,29]), and by Jordan G. Brankov and Valentin Zagrebnov in [31], dealing with the Husimi-Temperley model (and then deepened in e.g., [11,12,15,30]).
• Each binding site may or may not be occupied by a ligand: this allows us to code its state (empty versus full) by a Boolean variable. For the generic i th site, we will use an Ising spin σ i = ±1, where σ i = −1 represents an empty i th site, viceversa σ i = +1 means that the i th site is occupied. Clearly, if there are overall N binding sites, i ∈ (1, ..., N ).
• It is rather inconvenient, an ultimately unnecessary, to deal with the whole set σ i , ..., σ N if we are interested in the properties of large numbers of these variables (i.e., in the so called thermodynamic limit corresponding to N → ∞). If we want to distinguish between a fully empty state σ i = −1∀i ∈ (1, ..., N ) (ordered case), a fully occupied state σ i = +1∀i ∈ (1, ..., N ) (ordered case) or a completely random case where σ i = ±1 with equal probability (disordered case), it is convenient to introduce the order parameter for these variables as the magnetization M (this term stems from the original application of the Statistical Mechanics model in the context of magnetism) that reads as the arithmetic average of the spin state, namely There is a univocal relation between the magnetization in Physics and the saturation function in Chemistry, where, we recall, we denote with Y ∈ [0, 1] the fractional occupation of the binding sites. In fact, one has [7,8] Eq.(5) constitutes the first bridge between the Chemistry we aim to describe (via the saturation function Y ) and the Physics that we want to use (via the magnetization M ).
• All the binding sites interact with the ligand by the same strength. This is a standard assumption in Chemical Kinetics [8,50,65] and it means that the diffusion of the ligands is fast enough to ensure a homogeneous solution. The concentration of free ligands is mapped into a one-body contribution H 1 in the cost-function. This term encodes for the action of an external magnetic field in such a way that, if the field acting on the i th is positive, the spin will tend to align upwards (namely this direction is energetically favored), and vice versa. This homogenous mixing assumption translates into a homogeneous external field h, and the related contribution reads as Notice that h plays as a chemical potential and, consistently, it can be related to the substrate concentration as S 0 being the value of the ligand concentration at half saturation.
Eq. (7) constitutes the second bridge between the Chemistry we aim to describe (via the ligand concentration S) and the Physics that we want to use (via the magnetic field h).
• The binding sites can cooperate in a positive manner: this can be modelled by introducing a coupling between the σ variables. The simplest mathematical form is given by a two-body contribution H 2 in the cost function. This term encodes for the reciprocal interactions among binding sites and it reads as where J ≥ 0 is the interaction strength and the sum runs over all possible pairs; the normalization factor 1/N ensures the linear extensivity of the cost-function with respect to the system size. A positive value for J implies an imitative interaction among binding sites: configurations where spins tend to be aligned each others (namely where sites tend to be either all occupied or all unoccupied) are energetically more favoured and will therefore be more likely.
• Combining together the previous contributions we get the total Hamiltonian It is possible to introduce the free energy associated to such an Hamiltonian as where β is the inverse temperature in proper units. The free energy is a key observable because it corresponds to the difference between the internal energy U and the entropy S (at given temperature), i.e. A(β, J, h) = S(β, J, h) − βU (β, J, h). If we could obtain an explicit expression for A(β, J, h) in terms of the order parameter M , we could obtain an expression for the magnetization expected at equilibrium by imposing δ M A(β, J, h) = 0, in fact, this implies that we are simultaneously asking for the minimum energy and the maximum entropy.
Notice that, having stated the two bridges given by Eq.s (5) and (7), other mappings between the two fields (e.g., the relation between the coupling strength J and the Hill coefficient n H , see eq.(25) later on) emerge spontaneously as properties of the thermodynamic solutions of the problems.
B. Resolution of the problem: the mechanical action We want to find an explicit expression (in terms of M ) for the free energy defined in eq. (10). To this task let us rename −βJ = t and βh = x and let us think of these fictitious variables as a time and a space, respectively. Thus, we can write the free energy as where we also wrote i<j σ i σ j ∼ (1/2) i,j σ i σ j , which implies vanishing corrections in the thermodynamic limit. If we work out the spatial and temporal derivatives of the free energy (12) we obtain where the average · t,x for a generic observable O depending on the spin configuration is defined as and, posing t = −βJ and x = βh, the Boltzmann average for the original system (9) is recovered and this shall be simply denoted as · If we now introduce a potential V (t, x), defined as half the variance of the magnetization, i.e., we see that, by construction, the free energy of this model obeys the following Hamilton-Jacobi equation and therefore A(t, x) is also an action of Classical Mechanics. We can simplify the previous equation by noticing that, for large enough volumes, the magnetization is a self-averaging quantity [14,63], thus in the infinite volume limit the potential must vanish, that is, lim N →∞ V (t, x) = 0. Here, we are restricting to large volumes and we are therefore left with a Hamilton-Jacobi equation describing a free propagation; since the potential is zero, the Lagrangian L coupled to the motion is just the kinetic term that is the analogous of the classical formula L = 1 2 mv 2 where the mass m is set unitary (i.e., m = 1), and the role of the velocity v is played by the average magnetization M . Solving the Hamilton-Jacobi equation is then straightforward: the solution is formally written as The evaluation of the Cauchy condition A(t = 0, x = x 0 ) is trivial because, at t = 0, the coupling between variables disappears (see eq.(10)), while the integral of the Lagrangian over time reduces to the Lagrangian times time (as the potential is zero). Pasting these two contributions together we obtain Finally, noticing that the equation of motion is a Galilean trajectory as recasting the solution back in the original variables, i.e. t = −βJ and x = βh, we get the free energy associated to this general positive cooperative reaction: By extremizing A(β, J, h) with respect to M we get This result recovers the well known self-consistency equation for the order parameter of the Curie-Weiss model in Statistical Mechanics [14,63].

C. Chemical properties of the physical solution
The self-consistent equation in Eq. (22) is an input-output relation for a general system of binary elements, possibly positively interacting, under the influence of an external field: the input in the system is the external field h and the output is the magnetization M . We can rewrite Eq.(22) in a chemical jargon by using the bridges coded in the Eq.s (5) and (7) and fixing, for the sake of simplicity, S 0 = 1, that is Before proceeding, we check that if cooperation disappears (i.e., binding sites are reciprocally independent), the Michaelis-Menten scenario is recovered. Posing J = 0 in the equation above we get that is (apart a constant factor that can be re-introduced by taking S 0 = k, rather than 1) the Michaelis-Menten equation (see Eq. 2). One step forward, we now take into account the coupling J and relate it to the Hill coefficient n H . The latter is defined in Chemistry as the slope of Y (S) at half saturation (i.e., when Y = 1/2), and we can obtain its expression following this prescription by using Eq.(23), namely We note that as J → 0 we get, as expected, n H → 1: if there is no cooperation between binding sites, the Hill coefficient must be unitary; further, the stronger the coupling J and the (hyperbolically) larger the value of the Hill coefficient. In particular, for J → 1 the kinetics gets ultra-sensitive and discontinuities emerge. We remark that, with simpler statistical mechanics model as linear chains of spins, phase transitions are not allowed, hence ultra-sensitive behavior can not be taken into account: the present framework is the simplest non-trivial scheme where all these phenomena can be recovered at once (see Fig. 1 and [7] for more details on ultra-sensitive kinetics). Also, it is worth highlighting the full consistency between our treatment of ultra-sensitive kinetics and more standard ones as for instance reported in [48] (see eq.5.17 therein), where the expression for the Hill coefficient can be translated into our formulation as One we see that for M → ±1 the Hill coefficient diverges, which is the signature of an ultra-sensitive behavior: this is perfectly coherent with our approach where, in that limit, the input-output relation (see the hyperbolic tangent (23)) becomes a step function. However, as mentioned in the Introduction, this theory has its flaws, in Chemistry as well as in Mechanics. Regarding the former, the complex picture of yeast's enzymes evidenced by Koshland [24,25], where positive and negative cooperativity appear simultaneously (and with the anti-cooperativity effect getting more and more pronounced as the substrate concentration is raised), still escapes from this mathematical architecture. Further, from the mechanical point of view, two weird things happen: the velocity M is bounded by c = 1, while in Classical Mechanics the velocity may diverge; further, if we look at the Boltzmannfaktor in the free energy (Eq. 12), this reads as exp N −tM 2 /2 + xM and, recalling that the kinetic energy in this mechanical analogy reads as M 2 /2 (the mass is unitary, thus velocity and momentum coincide), we are allowed to interpret A(β, J, h) as a real action. From this perspective, the exponent can be thought of as the coupling between the stress-energy tensor and the metric tensor: a glance at the form of the Boltzmannfaktor reveals that the natural underlying metric is (−1, +1) rather than (+1, +1) as in classical Euclidean frames, or in other words, it is of the Minkowskian type. All these details point toward the generalization of the equivalence including special relativity. Plan of the next section is to follow the mechanical path and extend the classical kinetic energy including relativistic corrections and then to investigate its implications. We will see that in the broader, relativistic framework for chemical kinetics the deviations that Koshland explained adding an anti-cooperative interactions -beyond the cooperative onesat high ligand's doses are the chemical analogies of the deviation from classical mechanics at high velocities observed in special relativity.

A. Relativistic setting
The relativistic extension of the the Hamiltonian (9) is defined by an Hamiltonian of the form where M = 1 N N i σ i as usual. Next, we have to insert (26) into the free energy (10): where we already replaced t = −βJ and x = βh in order to work out their streaming, that read as where the Boltzmann averages . . . t,x are defined as (using the magnetization as a trial function) As before, the averages . . . t,x will be denoted by . . . whenever evaluated in the sense of thermodynamics (i.e. for t = −βJ and x = βh). By a direct calculation, it is straightforward to see that the expression (27) obeys the relativistic Hamilton-Jacobi equation where the symbol represents the D'Alambert operator and V (t, x) is the potential whose expression is chosen in order to make the equation valid by construction and, this time, it is automatically Lorentz invariant. If the functional A(t, x) is sufficiently smooth (that is, its derivatives are regular functions of t and x), in the thermodynamic limit, we have lim N →∞ V (t, x) = 0, hence in this high-volume limit we are left with which is the Klein-Gordon equation for a free relativistic particle with unitary mass in natural units (m 0 = 1). In relativistic mechanics, the stress energy tensor of this particle is defined as where v is the classical velocity of the particle, γ = 1/ √ 1 − v 2 and E = γm 0 = γ is the relativistic energy. In addition, the contravariant momentum is expressed through the action by the following equation Comparing (32) and (33), it is immediate to identify the magnetization as the relativistic dynamical variable: while the Lorentz factor is In the thermodynamic limit, the particle is free and its trajectories are the straight lines x = x 0 + vt. Since the relativistic Lagrangian L = −γ −1 is constant along these classical trajectories, the free energy can be computed as: Setting t = −Jβ and x = βh, we finally get an explicit expression for the free energy: Requiring that the free energy is extremal with respect to the magnetization (that from a thermodynamical perspective can be seen as the simultaneous effect of the minimum energy and the maximum entropy principles, and, from a mechanical perspective as the minimum action principle), the associated self-consistency equation becomes

B. The classical limit from a chemical perspective
Reading the self-consistency (38) in chemical terms, that is using the bridges Eq.s (5,7), we obtain We can now check whether, under suitable conditions, this broader theory recovers the classical limit. First, we notice that under the assumption of no interactions among binding sites (i.e., J = 0) and replacing h = 1 2 log(S/S 0 ), the Michaelis-Menten behaviour is recovered. This can be shown by rewriting Eq. 39 as where we also shifted S/S 0 → S for simplicity. For J = 0 the previous equation reduces to Y (S) = S/(1+S). Further, taking the classical limit, at the lowest order, we have the following expansions such that (38) reduces to (22), in the physical context, and to (23), in chemical context. Clearly, also the slope at Y = 1/2 is preserved hence, in the classical limit, we recover the expected expression for Hill coefficient (see Eq. 25), namely

C. Beyond the classical limit
To understand why we expect variations with respect to the Hill paradigm at relatively large values of the substrate concentration, we must check carefully the relativistic self-consistency (38). Let us assume we are working at nottoo high velocities (that is M < 1) and we can expand the argument inside the hyperbolic tangent, in particular, approximating 1/( The relativistic effects in chemical kinetics become transparent in this way: if we look at the field felt by the binding sites (i.e., the argument inside the hyperbolic tangent), we see that, beyond the standard Hill term βJ M (that positively pairs binding sites together) another term appears that, this time, negatively pairs binding sites together. Retaining this level of approximation, we could write an effective Hamiltonian to generate Eq. (43) that reads as hence, beyond the two-body positive coupling coded by the first term, another four-body negative coupling appears. The latter is responsible for the deviation from the classical paradigm and these deviations are in full agreement with the Koshland generalization toward the concept of mixed positive and negative cooperativity [24]. In particular, we can see at work the entire reasoning of Koshland who pointed out how, at large enough substrate concentration, the positivity of the reaction diminishes. In fact, for M ∼ 0 no relativistic effect can be noted. By increasing S (the input in the system), we get a growth in M (the output in the system): the latter raises in response of S and it is enhanced because of the two-body term in the effective Hamiltonian (44), the four-body term still being negligible. As S keeps on growing, M increases as well, up to a point where it reaches high enough values such that, from now on, also the four-body term inside the effective Hamiltonian (44) becomes relevant. At this point, a novel, anti-cooperative effect is naturally induced in the reaction and it yields to a reduction of the Hill coefficient. In the next analysis these qualitative remarks shall be addressed in more details. We focus on the definition of the Hill coefficient based on the Hill equation This equation accounts for the possibility that not all receptor sites are independent: here n H is the average number of interacting sites and the slope of the Hill plot. The latter is based on a linear transformation made by rearranging Eq. (45) and taking the logarithm: Thus, one plots log Y /(1 − Y ) versus log S, fits with a linear function and the resulting slope, calculated at the halfsaturation point, provides the Hill coefficient. As already underlined, the Michaelis-Menten theory corresponds to n H = 1 and any deviations from a slope of 1 tell us about deviation from that model. For the (classical and relativistic) models analyzed here (coded in the Hamiltonians (9) and (26)) we can estimate the slope n H directly from the self-consistency equations (23) and (39). Let us start with the classical model. We preliminary notice that Therefore, we just need to evaluate dY /d log(S) in Y = 1/2, which reads as Posing x = dY /d log(S)| Y =1/2 and noticing that S = 1 when Y = 1/2, we have By plugging this result in Eq. 47, we finally have One can see that when J = 0 the Hill coefficient is unitary as expected for non-cooperative systems, when J > 0 the coefficient is larger than 1, indicating that receptors are interacting, and when J < 0 the coefficient is smaller than 1, as expected for negative cooperativity.
Let us now move to the relativistic model. Again, we just need to evaluate dY /d log(S) in Y = 1/2, which, recalling (39), reads as Exploiting the fact that S = 1 when Y = 1/2, the previous expression simplifies as Thus, we can write Note that n class H /n rel H < 1, confirming that the relativistic correction weakens the emerging cooperativity.

D. Further robustness checks
As stressed above, for a fixed interaction coupling J, the relativistic model is expected to exhibit a lower cooperativity with respect to the classical model. In order to quantify this point we considered different quantifiers for cooperativity and we compared the outcomes for the relativistic and the classical models set at the same value of J. Let us start with the Koshland measure of cooperativity which is defined as the ratio [66]  where R 0.9 denotes the substrate concentration corresponding to a 90% saturation, while R 0.1 denotes the substrate concentration corresponding to a 10% saturation, that is, Y (R 0.9 ) = 0.9 and Y (R 0.1 ) = 0.1. In the non-cooperative case one has R 0.9 /R 0.1 = 81 and, accordingly, if the ratio is smaller than 81 (meaning that the saturation curve is relatively steep) one has positive cooperativity, while if the ratio is larger than 81 one has negative cooperativity. The advantage in using the index κ is that it can be easily measured, yet it ignores all information that can be derived from the shape of Y (S). In particular, this quantifier can be estimated starting from a Klotz plot (see e.g., Fig. 2, panel a) where the saturation function is shown versus the logarithm of the (free) ligand concentration; in the presence of positive cooperativity this plot yields to a characteristic sigmoidal curve. For the models analyzed here we can estimate R 0.9 /R 0.1 directly from the self-consistency equations (23)- (39). Starting from the classical model and, posing Y = 0.9 and Y = 0.1 we get, respectively Of course, when J = 0 we recover the value 81, when J > 0 we get R class < 81, and when J < 0 we get R class > 81. Repeating analogous calculations for the relativistic model we get and, with some algebra, log(S 0.9 ) = 2 atanh that is, Again, one can check that when J = 0 we recover the value 81, when J > 0 we get κ rel < 81, and when J < 0 we get κ rel > 81. Also, κ rel /κ class = e −16J/ √ 41+16J/5 > 1. This means that, even with this quantifier, when fixing the same coupling constant J, the emerging cooperativity is weaker for the relativistic model, as expected.
Next, let us consider the cooperativity quantifier derived from the Scatchard plot. We recall that this plot is built by showing the behavior of Y /S with respect to Y . In fact, according to the simplest scenario [67], at equilibrium, one can write where k is the proportionality constant between response and occupancy (i.e., it is the ratio between the dissociation and the association constants), and rearranging Eq. 62 we have The previous expression fits the equation of a line for Y /S versus Y , whose slope is −1/k. The advantages in using the Scatchard plot is that it is a very powerful tool for identifying deviations from the simple model, which, without deviations, is represented by a straight line. In particular, a concave-up curve may indicate the presence of negative cooperativity between binding sites, while a concave-down curve is indicative of positive cooperativity. Also, in the latter case, the maxima occurs at the fractional occupancy Y * which fulfills where σ provides another quantifier for cooperativity.
Starting from the classical model, we can build the function Y /S, by first getting S as a function of Y and can be obtained by inverting formula (23), namely By deriving Y /S with respect to Y we get which is identically equal to −1 when J = 0, monotonically decreasing with Y when J > 0 and monotonically increasing with Y when J < 0. The (possible) root therefore provides the extremal point, that is and, comparing with Eq. 64, we get We now repeat analogous calculations for the relativistic model. First, we get S as a function of Y , by inverting formula (39), namely By deriving Y /S with respect to Y we get d dY We can obtain an estimate of the value Y * corresponding to the maximum by recalling Y ≤ 1 and neglecting high-order terms. In this way we get and, comparing with Eq. 64, we get The three plots considered here (i.e., Klotz, Scatchard, and Hill), and the related estimates for the extent of cooperativity, are presented in Fig. 2. In particular, in panel d) we compare the cooperativity quantifiers for several values of J: as anticipated, in general, for a given value of J, the relativistic model gives rise to a weaker cooperativity.
We proceed our analysis by deepening the role of the coupling constant J in the binding curves related to the two models. In Fig. 3 we present the Klotz's plot (panel a), the Scatchard plot (panel b), and the Hill plot (panel c) for the relativistic and the classic models, comparing the outcomes for different values of J. As expected, the point corresponding to S = 1 and Y = 1/2 is a fixed point in each plot and, in general, the gap between the two models is enhanced when J is larger (i.e., when J is closer to 1). Also, when J is not too small, the Scatchard plot for the relativistic model displays a flex at small values of Y suggesting that, when the saturation is small, the cooperativity is not truly positive.
In the final part of this section we want to get deeper in the comparison between the classical and the relativistic models. To this aim, we solved numerically Eq. 39, for different values of S and of J, getting a set of data Y (S, J). We can think of this set of data as the result of a set of measurements where we collect the saturation value at a given substrate concentration. Now, assuming that in this experiment we have no hints about the underlying cooperative mechanisms, we may apply the formulas for the plain positive cooperativity and infer the value of J. More practically, we calculate numerically Y from the relativistic model for different values of S and of the coupling strength, referred to as J rel for clarity. Next, we manipulate the set of data Y (S, J rel ) by inverting the formula in Eq. 23: as the value of S is assumed to be known, we can derive the coupling strength, referred to as J class , expected within a classical framework. In this way, we can compare the original coupling constant J rel with the inferred one J class . We can translate this procedures in formulas as follows:  (39). From this data we inferred the expected classical coupling J class by inverting the self-consistent equation (23). We repeated the same operations for several values of S and J rel . In the left panel we show the inferred J class versus the fixed J rel : different colors represent different values of S and the identity function is also shown for reference (dashed, black curve). Notice that, in general J class < J rel and the inequality is enhanced as S grows. In the right panel we show a contour plot for the ratio J class /J rel versus h = log(S)/2 and J rel . Again, one can notice that, in general, J class /J rel < 1 and this inequality is enhanced for relatively large values of S.
In Fig. 4 (panel a) we plot J class versus J rel , for different values of S. Notice that the two parameters are related by a linear law, whose slope is smaller than 1 and decreases with S. This confirms that the relativistic model yields to a weak cooperativity. The negative contributions in the relativistic model get more effective when J rel and S are large, as further highlighted in Fig. 4 (panel b).

IV. CONCLUSIONS
The rewards in the overall bridge linking Chemical Kinetics and Analytical Mechanics are several, both theoretical and practical, as we briefly comment. The former lie in a deeper understanding of the mathematical scaffold for modelling real phenomena: it is far from trivial that the description of chemical/thermodynamical equilibrium is formally the same as the mechanical one. In particular, the self-consistency relation (38) that emerges from the thermodynamic principles (in fact, it stems from the requirement of simultaneous entropy maximization and energy minimization) also turns out to be, in the mechanical analogy, a direct consequence of the least action principles δA(t, x) = 0. This means that the stationary point corresponds to a light perturbation of the evolution of the system in the interval [0, t]. Explicitly, we shift infinitesimally M t,x → M t,x + δ M t,x , then 0 = δA(t, x) = ∂A(t, x) from which (38) is recovered (as usual by setting t = −Jβ and x = βh), since this holds for all variations δ M t,x . Even more exciting, still by the theoretical side, is the realization of the complexity of systems presenting mixed reaction (i.e., where both positive and negative cooperativity are simultaneously at work) and the possible applications in information processing, as we are going to discuss. First, let us clarify that in the Literature we speak of complex network or complex system with (mainly) two, rather distinct, meanings: in full generality, let us consider an Hamiltonian as H(σ, J) = i<j J ij σ i σ j and let us write the two-body coupling matrix as J ij = A ij W ij , where A is the adjacency matrix, accounting for the bare topology of the system (its entry A ij is 1 if there is a link connecting the related nodes (i, j), which are therefore allowed to interact each other, and it is zero otherwise) and W is the weight matrix, accounting for the sign and the magnitude of the links (i.e. the type of interactions among binding sites).
Dealing with A, networks where the topology is very heterogeneous (e.g., the distribution of the number of links stemming from a node has a power-law scaling) are called complex networks, as it is case for the Barabasi-Albert model [13]. Dealing with W , networks where the entries of the weight matrix are both positive and negative are termed complex systems, as the Sherrington-Kirkpatrick model [51] for the so called spin-glasses. Crucially, spin glasses spontaneously show very general information-processing skills and computational capabilities: for instance Hopfield neural networks [35] and restricted Boltzmann machines [1] -key tools in Artificial Intelligence (respectively in neural networks and machine learning)-are two types of spin-glasses and it is with this last definition of complexity that we now can read the information processing capabilities of the elementary reactions we studied. For a given macromolecule under consideration, we could paste each binding sites on a node and draw the links among nodes that are interacting: if two nodes are correlated (they show positive cooperativity), their relative interaction is positive, while if two nodes are anti-correlated (they show negative cooperativity), their relative interaction is negative. Dealing with mixed reactions we have to deal with spin glasses and we can thus assess how much information has been processed in a given reaction by evaluating the amount of information processed in its corresponding spin-glass representation, using our bridge. We have already started this investigation in [5][6][7]. Finally, from a practical perspective, in the classical limit (i.e., for simple reactions) we have an explicit expression that directly relates the Hill coefficient n H , which can be measured experimentally, to the interaction coupling J assumed in the model, namely n H = 1/(1 − J). This allows designing specific models and very simple validations (at least at this coarse-grained level) and gives a new computational perspective by which analyze already developed ones (see e.g. [26,33,37,40,47]. Then, regarding complex reactions, the puzzling scenario evidenced by Koshland, finally finds out a simple descriptive framework that, crucially, also recovers to the standard cooperative scenario in the proper limit: full coherence among various, apparently antithetic, results is obtained within a unique framework.
The author declare that there is no conflict of interest regarding the publication of this paper.