Feedback-induced self-oscillations in large interacting systems subjected to phase transitions

In this article it is shown that large systems with many interacting units endowing multiple phases display self-oscillations in the presence of linear feedback between the control and order parameters, where an Andronov–Hopf bifurcation takes over the phase transition. This is simply illustrated through the mean field Landau theory whose feedback dynamics turn out to be described by the Van der Pol equation and it is then validated for the fully connected Ising model following heat bath dynamics. Despite its simplicity, this theory accounts potentially for a rich range of phenomena: here it is applied to describe in a stylized way (i) excess demand-price cycles due to strong herding in a simple agent-based market model; (ii) congestion waves in queuing networks triggered by user feedback to delays in overloaded conditions; and (iii) metabolic network oscillations resulting from cell growth control in a bistable phenotypic landscape.


Introduction
From clocks to motors [1], many musical instruments [2,3] (including our throat [4]), pumping hearts [5], and firing neurons [6], many oscillatory systems are self-sustained [1]. Despite their many applications [7] to describe natural and engineering systems, a thorough physical understanding of self-oscillations is lacking [8], in particular their emergence as a collective behavior in large systems, like the ones studied by statistical mechanics. This in turn hampers thermodynamic analysis, an aspect that has puzzled scholars since the first seminal articles on the subject, where self-oscillations were provokingly regarded as a form of perpetual motion [9]. Oscillatory collective behaviors of many interacting units have been mainly investigated from the point of view of their synchronization where the interacting units are already assumed as linear oscillators [10,11] or self-oscillators [12,13]. In this work a different route is taken and a general criterion will be given for the onset of self-oscillations in large systems in a fully self-organized way, without postulating that the elementary units are oscillators. Rigorous criteria exist for the development of self-oscillations in dynamical systems [14]. From a more physical viewpoint, in particular in the context of electrical engineering, self-oscillations are triggered for a workload corresponding to the 'negative resistance' part of the current-voltage characteristic curve of active devices [1]. The key idea of this work is that systems with many interacting degrees of freedom that show phase coexistence, upon treating the equation of states on equal footing of characteristic curves, develop self-oscillations in the presence of feedback between the control and order parameters that try to force them on thermodynamically unstable branches. This is illustrated in a Gedankenexperiment in the next section where the feedback dynamics of the Landau mean field theory turn out to be described by the Van der Pol oscillator, a prediction that is successfully tested for the fully connected Ising model subject to heat bath dynamics. While in the Ising model the feedback is artificially introduced for illustrative purposes, in the following sections several different phenomena are considered in a stylized way in which the feedback is fully dynamically justified, namely (i) excess demand-price cycles due to strong herding in a simple agent-based market model; (ii) congestion waves in queuing networks triggered by user feedback to delays in overloaded conditions; (iii) metabolic network oscillations resulting from cell growth control in a bistable phenotypic landscape. In the conclusion results are summarized and several interesting outlooks that stem from this work are pointed out.

Self-oscillating Ising model
Consider the following simple gedankenexperiment: we have a magnet in equilibrium in a given external magnetic field h and thermal bath of temperature T < T c (where T c is the critical temperature for the para-ferromagnetic transition) and we measure its magnetization m(T, h). Can we control and set it to m = 0 by tuning the external magnetic field with simple linear negative feedback dh ∝ −mdt , at fixed temperature?
In the following it will be shown that this is not the case at least for the mean field case. In order to address the question formally, we consider the Landau expansion of the free energy density in the magnetization m [15] with parameters β = 1/T and h: Equilibrium relaxation dynamics approximately follow the phenomenological equation (linear response, where the time unit coincides with the typical one) [16] The extrema of the free energy are thus equilibrium states of the system: without an external magnetic field h = 0 we have m = 0 as the unique stable equilibrium point (node) as soon as β < 1, whereas for β > 1 this becomes unstable while the other two nodes appear at m = ± 3(β − 1) which are stable equilibrium points (spontaneous symmetry breaking). Let us now add the linear negative feedback between the external field and the magnetization with timescale τ , i.e. τḣ = −m. We have the dynamical systeṁ It is easy to see that for β > 1 this system has no stable steady states. Upon differentiating equation (3) and substituting the value of ḣ from equation (4) (Lienard transformation) we have the equivalent second-order system This is an instance of the Van der Pol equation [17], which is known to display self-oscillations for β > β c = 1 and the phase transition has been substituted by an Andronov-Hopf bifurcation. This prediction has been tested in the simplest microscopic setting, that is the fully connected Ising model. This is composed of N spin variables s i = ±1 interacting through pairwise homogeneous ferromagnetic couplings of strength J whose energy E in an external magnetic field h is ( M = i s i ) and whose equilibrium configurations are given by the Boltzmann-Gibbs distribution We will consider single spin flip heat bath dynamics, i.e. at each time step we choose one spin s i uniformly at random and set its new state, given the magnetization density m = M/N , with probability and update the external field with the prescribed feedback rule (∆t = 1/N) Analytical insights can be obtained by means of a Van Kampen system size expansion along the following lines. Upon considering the induced evolution for the magnetization M from  the heat bath dynamics of the single spins, this performs a random walk in the interval [−N, N] with step size −2, 0, 2 and the (normalized) rates

The master equatioṅ
can be expanded in the system size [18], i.e. upon considering a decomposition of M under a scaling hypothesis in terms of the auxiliary variables performing an expansion of the master equation in the parameter 1/ √ N and neglecting higher order terms, we have, upon considering on one hand a deterministic equation for φ on the other a linear Fokker-Planck equation for z where both depend on the external parameter h, which can be considered varying in time and subject to the negative feedback Monte Carlo simulations as well as analytical calculations obtained by the Van Kampen expansion clearly confirm the predictions, in particular • for β ∼ 1, the system shows approximately harmonic oscillations where the amplitude adjusts itself making the damping term negligible, i.e. m 2 ∼ β − 1 (figure 1).
• for β 1 it performs relaxation oscillations: approximately, in the phase space (m, h), dynamical trajectories follow the isocline curve (where dh dm → ∞) up to the tipping unstable points, where they follow straight horizontal lines corresponding to sudden jumps (figure 2). Furthermore, the Van Kampen approximation analytically highlights a typical trend for the variance that turns out to have its maxima in correspondence to the tipping points. A deterministic equation for this term can be worked out from the linear Fokker-Planck equation (16) by considering the equivalent Ito stochastic differential equation (where W(t) is a Wiener process) (18) and upon differentation from which we obtain the linear equation for the variance (σ 2 = z 2 ) This fluctuation term from the Van Kampen expansion has been compared to simulation results, finding a good agreement (see figure 3 below). The mean field approach (in particular equation (5)) predicts that the oscillatory time period scales as T ∝ √ τ for β β c = 1 and as T ∝ τ for β 1, in good agreement with numerical results (see figure 4).
It is worth noting that the onset of self-oscillations in the Ising model is not limited to infinite dimensional graphs. Monte Carlo simulations give evidence of self-oscillations in a 2D square lattice (figure 5) if feedback is provided between local external fields and single spins, i.e.
Although the feedback has been introduced artificially here for illustrative purposes, the analyzed dynamics could have interesting applications to the study of magnetic refrigeration cycles by magnetocaloric effect [19].

Excess demand-price cycle: herding in a simple agent-based model
The Ising model is a paradigm for the study of collective behaviors with applications extending beyond physics. In particular much research has been devoted in recent times to using concepts and methods from the statistical mechanics of the Ising model to analyze agent-based models (ABMs) for social and economic phenomena [20], in particular oscillatory behavior and synchronization transition in crisis waves by means of macroeconomic ABM [21,22]. In this section we will illustrate the emergence of business cycles [23] in a simple ABM as a consequence of the proposed theory. We will focus here on a simple model sharing similarity with minority games [24], where N agents shall choose one of two possible actions, i.e. the state of agent i is given by a spin variable s i = ±1. Upon identifying in a stylized way these two states as the propensity to sell or buy a given good, the magnetization M = N + − N − is called the excess demand and phenomenologically price dynamics shall follow (if more agents sell than buy, the price rises and vice versa) where the price has been called h in order to point out the analogies with the Ising model; this equation is thus equivalent to the feedback introduced in the previous section. Given a reference price that we will set for simplicity to zero, we will assume that agents update their state with a certain frequency ν stochastically with probability (bounded rationality) That takes into account the fact that agents have larger propensity to buy (sell) if they perceive that the price is below (above) the reference value. It can be easily seen that this simple model has a stable Gaussian fixed point. On the other hand, it is known that economical and social agents can often make decisions based on the imitation of neighbors in their social network. The cumulative effect of this imitation (herding) is known to account for instabilities in the underlying dynamics [25]. We will explore this issue of network embeddedness and herding by recasting the aforementioned model in terms of the co-evolving network models studied in [26], where social links are created and destroyed among the agents, in turn reinforcing their beliefs, specifically: • Links are created with rate η (per agent) among agents sharing the same state.
• Links are destroyed with rate 1.
We will consider a regime of strong herding (that is the zero temperature limit also considered in [26]), i.e.
• Agents update their state upon looking at the price only if they are isolated.
At fixed h the equation of the state of the system can be worked out approximately with the methods of [26]. The approximate population dynamics equations for the density of agents with k connections in state σ read: If we consider the generating function G σ (s) = k n k,σ s k we have in the steady state G ss σ (s) = n 0,σ e xσs and the self-consistent equations ν + n 0, x σ /η = n 0,σ e xσ (28) x If we parametrize ν − /ν + = e h and consider m = x+−x− η we have finally the equation of the state that shows a second-order critical point at (η c = 2, h c = 0). Upon considering phenomenologically the relaxation and considering the price feedback dynamics (with timescale τ ) we finally have the systeṁ This prediction is tested against Monte Carlo simulations where a good agreement is found (see figure 6). Interestingly, in this case the whole network of interactions is oscillating in time, assuming the form of a Poissonian tree-like random graph (Erdos-Renyi) whose average degree oscillates. This can be seen as an instance of a self-oscillating temporal network [27].

Congestion waves in queuing networks
In this section we will test the theory in the context of out-of-equilibrium systems. This is a very rich topic, difficult to study given the lack of general established variational principles, like maximum entropy, leading to the Boltzmann-Gibbs measure characterizing their equilibrium counterpart. Perhaps one of the simplest yet general class of processes displaying phase transitions in this context are the so-called zero-range processes [28]. Broadly speaking, these are models of particles hopping among nodes in a graph whose hopping rates depend only on the number of particles present on the departure and arrival nodes. This very short range interaction permits factorization of the steady-state probability that is amenable for analytical calculations [28]. Among them, we have queuing networks [29], which are known to be subject to a congestion phase transition [30,31]. In these systems, used to model communication networks in engineering studies, packets of information are injected into the network by users for processing purposes and stored in queues at the nodes. If the load that the network experiences overcomes its processing capabilities, queues will start to grow congesting the system. In a stylized way, upon calling n the number of packets in the network, p the packet injection rate ('load') and considering a function f (n) that encodes for the network processing rate, we have the average rate equation: The function f (n) depends on the chosen model, in general f (0) = 0, f 0, lim n→∞ f (n) = µ < ∞. If p > µ we will have a steady growing solution ṅ p − µ, while stationary solutions are given by n s = f −1 ( p) and are stable if f (n s ) > 0. Generally unstable fixed points come with coexisting stable ones in this setting. Traditional work on queuing network theory focused on the stationary regime while recent observations of traffic in large networks spurred an increasing interest in the congested regime. However, traffic data seem to show rather more complex phenomena with respect to simple steady growing queues, including waves and intermittent behaviors [32], and it has been pointed out that one of the key ingredients accounting for the latter relies on the level of user feedback to overloaded conditions [33]. In the most simple setting this is naturally implemented within our scheme by introducing linear feedback dynamics for the load p: users have typical reaction times τ , they strive for a desired load p d , and they reduce the load if queues are longer (linear negative feedback of strength c): In the presence of the feedback the dynamics could present no stable steady states but selfoscillations and it will be shown that this is the case in the presence of a simple rule of traffic control. Previous studies have shown that traffic control, mimicking known Internet protocols like TCP/IP, can enhance the processing capabilities of the system, increasing the free-flow region but at the price of introducing non-linearities that trigger congestion transition in a discontinuous way, with hysteresis and coexistence [34,35] that are reflected in a non-monotonous f (n). In correspondence of these points we expect the linear feedback to induce selfoscillations that shall be thus specific to cases in which traffic is controlled.
These predictions have been successfully tested on a Jackson or open queuing network, consisting of N nodes such that: • each node i is endowed with a first-in-first-out queue with unlimited waiting places (it can be arbitrarily long). • The delivery of a packet from the front of i follows a Poisson process with a certain frequency t i (service rates), and -the packet exits the network with some probability µ i , or -it goes on the 'back' of another queue j with probability q ij .
• Packets are injected in each queue i from external sources by a Poisson stream with intensity p i .
We considered random walk routing q ij = 1/k j where k j is the degree of the receiving node j, and completely homogeneous conditions p i = p, µ i = µ, t i = t = 1. Traffic control has been included with the following simple rule [34]: • The receiving node j starts to reject particles with probability η once its queue is longer than n * Results shown in figure 7 are obtained by simulations on an Erdos-Renyi random graph (average degree z = 5, size N = 10 3 ), desired load and absorbing rate p d = µ = 0.2, user feedback strength c = 2 · 10 −6 , and traffic control parameters η = 0.75, n * = 10. In the presence of feedback dynamics for p, we would expect now that upon loading the system by crossing f (n) in its decreasing part, this enforces self-oscillations and this is verified as we show in figure 8. It is worth noticing that the self-intersection of the trajectory in the plane (n, p) excludes by means of the Cauchy theorem a simple 2D mean field picture.

Autonomous metabolic network oscillations
In a biological context classical models in statistical physics are increasingly used to perform inference and analyze data: examples range from random walks and biopolymers [36], Ising models and neural networks [37] to continuous spin models and flocking birds [38] to cite just a few. Standard flux balance analysis (FBA) approaches to model cell metabolism [39] share a formal analogy with the Gardner problem in statistical mechanics [40] and it has been recently shown that maximum entropy inference schemes [41] outperform FBA in modeling flux data of the catabolic core of E. coli [42]. The metabolism is the network of enzymatic reactions that sustains the free energy needs of the cell, strongly constrained by physico-chemical laws that in turn provide for suitable modeling. Metabolic dynamics give well-known examples of nonlinear self-oscillators in particular glycolytic oscillations [43], experimentally tested in living cells [44]. With regard to whole cell metabolism, recent findings in yeast indeed show intrinsic whole single-cell metabolic oscillations, autonomous from the cell cycle and potentially able to drive it [45]. In this section we will explore, within the proposed theory, the possibility that such oscillations are in general due to the effect of feedback in the presence of a bistable phenotypic landscape. Feedback mechanisms are needed in order to maintain cell size homeostasis [46] and a bistable landscape has been observed in E. coli [47,48], where it is at the core of persistence phenomena [49]. Furthermore it has been recently pointed out theoretically within the framework of constraint-based models that second-order moment constraints on the growth rate in general enable bistability [50]. These are stationary models of large chemical networks including in a realistic way the stoichiometry of known pathways and a phenomenological biomass growth reaction that is a function of the enzymatic fluxes λ = λ(v). For a chemical reaction network in which M metabolites participate in N reactions with the stoichiometry encoded in a matrix S = {S µr }, the concentrations c µ change in time according to mass-balance equationṡ where v i is the flux of the reaction i (that is in general a function of the concentration levels v i (c)). The steady state implies S · v = 0. In constraint-based modeling, apart from mass balance constraints, fluxes are bounded in certain ranges v r ∈ [v min r , v max r ] that take into account thermodynamic irreversibility, kinetic limits, and physiological constraints. The set of constraints defines a convex polytope P in the space of reaction fluxes. We seek the states fixing the first two moments of the growth rate, given by the Boltzmann distributions: They have been sampled by means of a hit-and-run Monte Carlo Markov chain with ellipsoidal rounding [51] on the model of the catabolic core of E. coli from the genome scale reconstruction [52], including glycolysis, the pentose phosphate pathway, the Krebbs cycle, oxidative phosphorylation, nitrogen catabolism, and the biomass growth reaction, for a total of N = 95 reactions and M = 72 compounds, simulated in a glucose-limited aerobic minimal medium with a maximum glucose uptake of u g,max = 10 mmol/gdwh. Upon counting the number of feasible states leading to the same growth rate by uniformly sampling P, the marginal growth rate distribution can be recast in terms of the rate function F(λ) (where we posed λ max = 1, the maximum growth rate in the model obtainable by linear programming, and we get a simplex-like entropic term with a 20 for the carbon catabolic core of E. coli) It should be noted that such maximum entropy distribution can be the steady state of suited population dynamics, which for the case γ = 0 has been shown to be the logistic [53]. If we consider relaxation dynamics in the linear response regime and linear control through β, looking for a desired λ s , we have the dynamical systeṁ that upon Lienard transformation is mapped into the second-order system A sufficient condition to get self-oscillations is γ > a 2(1−λs) 2 . This is confirmed by simulations on a model of the catabolic core of E. coli (figure 8), where growth rate oscillations entrain correlated metabolic fluxes. Here, the activity of the enzyme cytochrome b oxidase and the nitrogen and phosphate uptakes are shown.

Conclusions
In this work it has been shown that self-oscillations arise in trying to control large systems that have been subjected to phase transitions. The oscillatory collective behavior emergent in this way does not depend on postulating oscillatory units, but it is fully self-organized. The key idea is that large systems endowing multiple phases develop self-oscillations in the presence of feedback that tries to force them on thermodynamically unstable branches, analogously to active electrical devices controlled with a workload corresponding to 'negative resistance' parts of their characteristic curve. It has been shown that the feedback maps the Landau mean field theory into the Van der Pol oscillator and such behavior has been confirmed for the fully connected Ising model subject to heat bath dynamics. Further studies are needed to test the theory on finite dimensions as well as experimentally. Preliminary numerical results on a 2D square lattice seem to show that self-oscillations are triggered by local feedback, while a global one triggers separation of magnetic domains. The theory is general and here it has been applied to describe in a stylized way (i) excess demand-price cycles due to strong herding effects in a simple agent-based market model; (ii) congestion waves in queuing networks triggered by user feedback to delays in overloaded conditions; (iii) metabolic network oscillations resulting from cell growth control in a bistable phenotypic landscape. All of these would deserve further work on their own, in particular suited analysis to test them against data by means of promising independent component analysis-based methods [54,55]. Apart from Andronov-Hopf, higher order bifurcations could be eventually connected to higher order terms of the self-interacting scalar field in the Landau approach and in turn connected to minimal lattice models whose interacting variables are ruled by less simple symmetry groups: tricritical Ising, Potts, etc. This would be a very exciting research program that would deserve further investigations. The theory could open the way to explore the thermodynamics and statistical physics of self-oscillations, and recent tools developed within stochastic thermodynamics could play a key role in this respect [56,57]. Finally, within the framework of control theory self-oscillations can be seen as an unwanted negative side effect and this theory could help shed light on the origin of this problem while controlling large complex networks [58].