Negative differential response in chemical reactions

Reaction currents in chemical networks usually increase when increasing their driving affinities. But far from equilibrium the opposite can also happen. We find that such negative differential response (NDR) occurs in reaction schemes of major biological relevance, namely, substrate inhibition and autocatalysis. We do so by deriving the full counting statistics of two minimal representative models using large deviation methods. We argue that NDR implies the existence of optimal affinities that maximize the robustness against environmental and intrinsic noise at intermediate values of dissipation. An analogous behavior is found in dissipative self-assembly, for which we identify the optimal working conditions set by NDR.


I. INTRODUCTION
Systems in contact with multiple (e.g. chemical, thermal) reservoirs fall out of equilibrium, in a state characterized by sustained mean currents (e.g. of matter, energy) [1]. These are controlled by affinities, the thermodynamic forces which measure the difference between the equilibria that distinct reservoirs try to impose on the system [2]. A perturbation in an affinity A-be it the deliberate manipulation of an experimenter or some environmental noise affecting the reservoirs-produces a small variation in a current J , quantified by the differential response function R = d J dA [3]. Close to equilibrium, such response is severely constrained [4]. Since currents are proportional to affinities, J = RA, the response R must be positive to ensure positivity of the entropy production Σ = J A = A 2 R 0 [5]. Far from equilibrium, instead, J need not be linear in A thus making R not only dependent on the entropy production. Kinetic aspects become relevant [6], thus opening the way to regimes of negative differential response (NDR) [7]. This counterintuitive, yet common phenomenon has been found in a wealth of physical systems after its first discovery in low-temperature semiconductors [8]. Examples are particles in crowded and glassy environments [9][10][11][12][13], tracers in external flows [14,15], hopping processes in disordered media [16,17], molecular motors [18,19], polymer electrophoresis in gels [20], quantum spin chains [21], graphene and thermal transistors [22,23]. The shared feature underlying all these systems is a trapping mechanism arising by (e.g. energetic, geometric, topological) constraints on the system states [24].
Here, we show that NDR plays a key role in open chemical reactions networks [25][26][27]. We show for three paradigmatic models-substrate inhibition, autocatalysis and dissipative self-assembly-how it appears in the average macroscopic behavior as well as in the stochastic regime. While the first two are well described core reaction schemes in living organisms [28,29], the latter is currently drawing the attention of chemists [30,31]. Within the scope of these examples we discuss the role of NDR with respect to environmental and intrinsic noise [32][33][34][35]. We first show that the region of marginal stability, i.e.
where R 0, ensures robustness against external perturbations (in the affinity) at moderate values of dissipation. We then argue that those systems affected by NDR that are not poised in the region of marginal stability, behave so in order to minimize the dispersion of the current. Such precision is found to be achieved at moderate values of dissipation, yet again. Hence, our findings show that the performance of life-supporting processes does not always increase at larger dissipation rates [36,37]. This rich behavior brought about by far from equilibrium conditions cannot be anticipated solely on the basis of general results, such as the recently derived thermodynamic uncertainty relations [38][39][40][41][42][43]. To unveil these properties, one needs to solve for the full counting statistics through large deviation methods [44]. Finally, since both robustness and precision are desirable in artificial applications of dissipative self-assembly, we identify the optimal affinity set by NDR using stochastic simulations.

II. THEORY
Because cells work at relative high, yet finite number of molecules, reaction currents fluctuate around their macroscopic average values. We assume the reactions to take place in a large well-mixed volume of size V , so that concentrations obey mass-action kinetics. The randomness of the single reaction events is described by the chemical master equation [45], that evolves the probability P t (c) of finding the concentration c σ of the dynamical species σ. Here, σ labels the dynamical species, while species whose concentration are fixed are labelled by σ . with which there occurs the reaction ρ involving ν σ,ρ (resp. ν σ ,ρ ) molecules of dynamical (resp. fixed) species σ (resp. σ ). The stoichiometric coefficient S σ,ρ = ν σ,−ρ − ν σ,+ρ then gives the net number variation of species σ per reaction ρ. To analyze the system response it is sufficient to focus on a reduced description based on the instantaneous number of reactions ρ per unit time, C ρ , and the typical rate to leave a chemical state c through reaction ρ, W ρ (c) (see Appendix). The complete statistics of their time-averaged value [47],(.) := 1 T T 0 dt(.), is encoded in the scaled cumulant generating function that gives upon differentiation all the covariances, e.g.
The averages . . . are performed along stochastic realizations with the path weight obtained from (1), that contains the auxiliary variable p accounting for random variations in particle number [48]. A standard technique to calculate (3) consists in absorbing the exponential counting factor of (3) into (4), changing H into the 'tilted' generator In view of the extensivity in T and V of the observables, averages performed with P q,λ are entirely dominated by the overwhelmingly more probable trajectory that maximizes (5). This observation allows us to calculate the scaled cumulant generating function as where c * and p * are solution of the steady-state Hamiltonian equations ∂ c H q,λ = 0 = ∂ p H q,λ . Currents can then be obtained as the net fluxes between forward and backward reactions,J ρ := (C +ρ −C −ρ ). The nonequilibrium origin of NDR emerges clearly from the stochastic setup. Indeed, the differential response of a generic current J ρ , can be obtained by expanding the generator H, and thus the path weight (4), to leading order in a small variation of the fixed concentration c σ . In general, it reads (see Appendix) whereρ are the reactions whose rates Wρ depend explicitly on the perturbed species σ . In (8) the current J ρ correlates with three distinct observables: the reaction currentJρ; the reaction trafficFρ := 1 T T 0 dt(Cρ + C −ρ ), i.e. the total unsigned number of ±ρ reactions; the reaction ratesWρ. Differently from currents, traffic and reactions rates do not have a definite thermodynamic character, their values being affected by kinetic factors. In the following we will focus on perturbations that alters only the af(see Appendix)finity A that drives J ρ . If such a perturbation happens at equilibrium, (8) reduces to the fluctuation-dissipation relation where only the entropic term R ∝ J 2 ρ cc appears (see Appendix). Out of equilibrium, instead, (8) shows that NDR arises when the currentJ ρ becomes sufficiently anticorrelated with either −Wρ orFρ. These two scenarios find their counterparts among physical systems undergoing mechanical trapping induced, respectively, by geometric constraints-a colloidal particle pulled through an array of obstacles [7] -and by many-body clustering-the same pulling experiment performed in a high-density medium [9].

III. SUBSTRATE INHIBITION
Substrate inhibition is estimated to occur in 20% of known enzymes [49]. In its simplest form [see fig. 1 (a)], it happens when up to two substrate molecules S can bind the active site of one enzyme E giving an inert species ESS. The binding of a single substrate molecule results in the formation of the active complex ES decaying into the product P, as in the usual Michaelis-Menten scheme. The latter pathway is responsible for the production of P from S at a concentration rate J 1 , that is the chemical current of biological interest. The former instead represents the competing process [28,49]. It takes upor traps, within the mechanical analogy-substrate into ESS thus decreasing the rate of production of P for large [S] ( fig. 1). Indeed, with [S] [P] kept constant by particle reservoirs to mimic physiological conditions and fixing the reaction affinity A = log k1k2 [S] k−1k−2[P ] , the stationary current takes the non-monotonic form [28,49] Here [E] tot is the total concentration of enzyme and . The kinetics of the usual Michaelis-Menten scheme is retrieved setting k 3 = 0 ( fig. 1).
The first two scaled cumulants of the time-averaged currentJ 1 show the existence of a marginal affinity A * that marks the transition to a NDR regime, i.e. R < 0 for A > A * , where fluctuations VarJ 1 := J 2 1 cc peak. In the present model of substrate inhibition, − J 1 (W 1 +W 3 ) cc is the leading negative contribution in (8) for A A * ( fig.  1), confirming that ESS is a trapping state.
The existence of NDR has some crucial consequences. First, since R(A * ) = 0, J 1 varies little upon sizable variations of substrate concentration around [S](A * ). Second, since J 1 is not an injective function of A, a target mean current-e.g. required for optimal physiological functioning-is attainable at two different affinities A min and A max . These two facts may constitute a crucial advantage to control environmental and intrinsic noise in biochemical systems.
In the first case, the system can reach a homeostatic state characterized by a relative stable output J 1 despite variations in the environmental conditions, i.e. the substrate concentration [S]. Importantly, a similar stable regime would be achieved only at larger affinities in the absence of NDR, i.e. for the standard Michaelis-Menten kinetics (cf. fig. 1). A representative example is the synthesis in neurons of dopamine (P) from tyrosine (S) mediated by the enzyme tyrosine hydroxylase (E) [50]. The tyrosine concentration in humans varies in response to meals on a timescale τ S ∼ 10 3 s, and typically ranges from 100 µM to 120 µM. Since the dynamics (1) for the substrate inhibition scheme in fig. 1 has a unique steady state, its typical relaxation timescale is well estimated by the inverse of the smallest (pseudo first-order) reaction . This interval is placed very close to A * 20.8, hence resulting in a current relative variation smaller than 3%.
In the second case, the system can increase the (scaled) signal-to-noise ratio SN R := J 1 / VarJ 1 selecting the optimal affinity among A min and A max . Consider, for example, the synthesis of serotonin (P) out of tryptophan (S) catalyzed by tryptophan hydoxylase (E) in human cells [51,52]. For different values of the parameters compatible with physiological conditions, we found that SN R is always smaller at A min , i.e. higher precision is achieved at A < A * (fig. 2). As a conse- yields a range of affinities 19.6 A 19.9 which is close to optimal in order to maximize SN R. Thanks to stochastic uncertainty relations [38][39][40][41][42], SN R can be bounded by dissipation, SN R Σ/2, and by the system's dynamical activity, The entropic bound means that a more precise current may be obtained at larger affinity, and thus dissipation. Nevertheless, such condition need not be realized in practice, especially because the bound becomes looser as A increases, as is the case for serotonin synthesis.

IV. AUTOCATALYSIS
Autocatalysis represents the second scenario in which NDR can arise, whose simplest possible scheme is depicted in fig. 3 (a). Having one dynamical concentration, two reactions (required to have a maximum current), and two fixed concentrations [S] and [P] (needed to set the system away from equilibrium), this is the minimal chemical scheme displaying NDR. An outstanding example falling into the autocatalytic paradigm is DNA replication [29]: two double stranded molecules are produced by one such molecule (X) and nucleobases (S), and eventually undergo a conformational change, e.g. into the double helix structure (P). Several other biological processes can be similarly described at a coarse-grained level as autocatalytic reactions, e.g. formation of micelles from amphiphiles [53,54], ATP net production in glycolysis [55], and conversion of prion proteins into the infectious form [56]. Here, we regard the autocatalytic scheme as a model for phosphorylation of protein kinase (X) coupled to a larger association/dissociation cycle via the con-version into the complex P [29,57,58]. For the chosen physiological values of the parameters, the coupled cycle is known to display circadian rhythmicity [58]. Taking [P] as time-independent, we highlight the role played by NDR in triggering chemical oscillations, a topic of major relavance which may even have a role in our understanding of the origin of life [59]. We consider the degradation of S into P as the current of interest. Its macroscopic value determined by the rate equation,  (8) shows that NDR is induced by a competition between forward and backward flows (due to the nonlinearity of the autocatalytic step), rather than by the presence of a trapping state. Also, the qualitative behavior of the VarJ 1 is different, with the minimum (rather then the maximum) occurring near A * . Despite that, for a wide range of parameters compatible with physiological conditions we observe that SN R is larger for A < A * , i.e. at smaller dissipation. Hence, as already discussed for substrate inhibition, autocatalysis can be run at low affinities to reduce the current dispersion or around the region of null response to minimize variations in the output current.

V. DISSIPATIVE SELF-ASSEMBLY
As a final example, we analyze dissipative selfassembly, a paradigm of out-of-equilibrium synthesis ex- tensively exploited by biological systems: prominent examples being the formation of microtubules out of tubulin dimers fueld by guanosine 5'-triphosphate (GTP) [60,61] and the ATP-driven self-assembly of actin filaments [62]. It has been also probed in experiments such as the controlled gelation of dibenzoyl-L-cysteine to form nanofibers [63] and the chemically fueled transient selfassembly of fibrous hydrogel materials [64]. A simple, yet insightful model is sketched in fig. 4 (a), which has been proposed as a minimal general scheme for genuine nonequilibrium self-assembly [31]. The direct aggregation of two monomers (M) to form the assembled state (A 2 )-which would be highly disfavored at equilibriumis boosted by coupling the process with the burning reaction of some fuel (F) converted into waste (W). This fueling mechanism opens side pathways involving the activated species M * , which easily aggregates into A * 2 . To give an example, supposing M to not aggregate because of unfavorable electrostatic interactions, then F (W) may be a high (low) energy methylating agent able to convert negatively charged monomers M into their neutral form M * . By properly fixing the concentrations of F and W, a nonequilibrium stationary state rich in the target species A 2 can be achieved. At odds with conventional equilibrium self-assembly, the efficacy of this synthetic procedure is not determined by the relative thermodynamic stabilities of the components, but rather by the sustained dissipation and kinetic aspects [30,65,66].
By design, the system attains large concentrations of A 2 depleting the monomer concentration [M] [31]. Therefore, the current of reaction ρ = 4-which is half the current from F to W-is almost unidirectional, especially far from equilibrium: Because of the proportionality relation (11), the existence of NDR affecting J 4 sets an upper bound on the maximal [A 2 ] achievable by the process. Being unable to calculate (3) for this model, we performed stochastic simulations based on the Gillespie algorithm [67]. We measured the mean current and its variance, as well as its response, for different values of the affinity A = log  (8) is not only conceptually revealing, but also of practical relevance for calculating responses without actually applying perturbations. Despite their proportionality in average, the currentJ 4 and the concentration [A 2 ] were found to possess different fluctuations. It implies that the signal-to-noise ratio [ is not bounded by dissipation, hence does not decrease at large A due to NDR. Indeed, [A 2 ] / Var[A 2 ] is close to its maximum at the optimal affinity A * . This is important for the scalability of artificial syntesis to microscopic volumes.

VI. DISCUSSION
In conclusion, we have shown that NDR is a widespread phenomenon in chemistry with major consequences on the efficacy of biological and artificial processes. In substrate inhibition, NDR allows a system to reach homeostasis at lower dissipation than in the Michaelis-Menten kinetics, keeping the signal-to-noise ratio unaltered. For systems that do not need to maintain a stable current, higher precision to sustain a given mean current can be reached at low affinity, i.e. dissipation.
Since the analogous behavior was found in both biochemical schemes, despite the difference in the qualitative behavior of the current fluctuations, the idea that life efficiency always increases with the dissipation rate is called into question [68]. Still, it is worth noticing that whenever these chemical schemes are used as effective models that coarse-grain some nonequilibrium reactions, the dissipation Σ is always smaller than the total entropy production rate of the original process [69]. Instead, if only equilibrated subprocesses are lumped or discarded, a complete thermodynamic description of the original process exists [70]. It identifies A with the chemical potential difference of the fixed species (respectively, P and S, F and W) and Σ with the entropy production rate [27]. Remarkably-given the scarcity of solvable models away from equilibrium-these results, obtained in the large-size limit, are exact. They show that the general bounds offered by the uncertainty relations have little predictive power for parameters that are biologically relevant.
Lastly, we have shown with stochastic simulations that NDR limits the efficacy of dissipative self-assembly: the ideal affinity that maximizes the output mean concentration also yields a nearly optimal signal-to-noise ratio. Altogether, we have achieved a fundamental analysis of NDR of reaction currents. It pinpoints the relation between robustness, precision and dissipation in biochemistry, and allows the optimization of performance and scalability in nonequilibrium synthesis.

VII. ACKNOWLEDGEMENTS
We thank Arthur Watchel for help with the numerical simulations, which were carried out using the HPC facilities of the University of Luxembourg [71]see https: //hpc.uni.lu.

Appendix A: Stochastic dynamics of chemical reaction networks
Consider a well-mixed volume V occupied by dilute reacting chemical species X σ , labelled by the index σ ∈ {1, . . . , M }, following mass action kinetics. The population number of the dynamical species n = (n 1 , . . . , n M ) varies in time because of the random reactions, while the concentration of the externally controlled species, c σ := [X σ ] with σ = {M + 1, . . . N }, is kept constant. A single reactive event, occurring thorough the reaction ρ ∈ {±1, . . . , ±M}, involves ν σ,ρ molecules σ and changes the population of species σ as n σ −→ n σ + S σ,ρ , with S σ,ρ := ν σ,−ρ − ν σ,+ρ the stoichiometric coefficient. For compactness, we will denote S ρ the vector of the stoichiometric coefficients corresponding to reaction ρ. The reactions happen with a probability rate with k ρ being the rate constant. The stochastic dynamics can be described by the chemical master equation that prescribes the time evolution of the probability P t (n) of the chemical populations in terms of the action of the operator H (V ) (n, ∂ n ). The solution of (A2) can be used to study only the statistics of state-like observables, i.e. functions of the instantaneous population n. In order to obtain the statistics of transition-like observables it is convenient to resort to a path integral representation of the probability of full stochastic trajectories. computed by functional integration of a path probability with 'tilted' generator The superscript (V, T ) stands for the dependence on a finite system volume V and trajectory duration T . Later, we will omit the superscripts V and T to indicate the large V and T limit of the various functions. All cumulants, such as the mean C (V,T ) ρ = ∂ qρ g (V,T ) (q, λ)| q,λ=0 and the connected correlations, e.g, −ρ ), follows from (A5) as well.

Macroscopic limit: the rate equations
In the thermodynamic limit, given by n σ → ∞, V → ∞ and finite concentrations c σ := n σ /V = [X σ ], the probability P t (n) becomes sharply peaked around its maximum. Thereby, one obtains the chemical rate equations multiplying (A2) by n and averaging,ċ with the average reaction current given by the (scaled) large-size limit of (A9): The same result can be obtained by the path integral formalism. The statistical weight in (A4) peaks as exp V T 0 dt[−ċ(t) · p(t) + H(c(t), p(t))] around those paths that maximize the time integral, i.e.,those that satisfy the Hamilton equationsṗ with H(c(t), p(t)) being now function of the macroscopic rates (A9). The rate equations (A8) are regain by looking for the noise-less trajectories: (A11)

Macroscopic limit: scaled cumulant generating function
Fluctuations in the thermodynamics limit are captured by the scaled cumulant generating function g(q, λ) := lim The average (A12) is dominated by the trajectories the maximize the statistical weight, namely, by the solutions of the Hamilton equationsṗ Since here we are interested only in systems with a single stable stationary state we focus on the unique fixed point of (A13)-we thus assume the absence of multiple stable fixed points and time-dependent attractors for (A13), which excludes the emergence in the thermodynamic limit of ergodicity breaking and limit-cycles, respectively. Namely, we seek the vectors c * and p * solution of To avoid clutter we do not explicitly write the parametric dependence of c * and p * on the counting fields q ρ and λ ρ . If there exist vectors λ such that λ · S ρ = 0 ∀ρ, then the dynamics (A2) conserves the concentrations λ · c. Therefore, (A13) needs to be supplemented by the constraints M σ=1 c σ λ σ = const, Using the solution c * and p * to evaluate (A12) yields the scaled cumulant generating function By virtue of the above assumptions, (A16) is a smooth function of q ρ and λ ρ .
Concerning the numerical values of the rate constants, for both the examples in the main text -i.e., tyrosine hydroxylase (TH) and tryptophan hydroxylase (TPH) -we relied on experimentally available K M = k2+k−1 k1 and K i = k−3 k3 [49]. Within these constrains, we chose realistic k ρ 's based on literature typical values [72]. In particular, we have set k 3 < k 1 , thus considering negative cooperativity between molecules S upon their binding to the enzyme E. k −3 has been kept small to make the "trapping effect" well highlighted, while k −1 and k 2 have been chosen in order to make reaction 2 rate limiting. k −2 is usually neglected in kinetic models, but here it guarantees thermodynamic consistency. Since the two enzyme considered have same K M and different K i , we have opted to keep differences minimal. Accordingly to the above argumentation, we checked the robustness of the qualitative features shown by the model under different choices. Plots in the main text were obtained with the following parameters: