Controlling a complex system near its critical point via temporal correlations

TL;DR

This paper investigates how complex systems near critical points can be controlled by leveraging temporal correlations, demonstrating a feedback mechanism that self-tunes systems like the Ising model, flocking, and neuronal networks to criticality.

Contribution

It introduces a novel feedback control method based on autocorrelation functions to steer diverse complex systems toward their critical points.

Findings

Feedback of autocorrelation functions can shift systems to criticality

Universal properties enable broad applicability of the control method

Validated on models: 2D Ising, 3D Vicsek, neuronal networks

Abstract

A wide variety of complex systems exhibit large fluctuations both in space and time that often can be attributed to the presence of some kind of critical phenomena. Under such critical scenario it is well known that the properties of the correlation functions in space and time are two sides of the same coin. Here we test wether systems exhibiting a phase transition could self-tune to its critical point taking advantage of such correlation properties. We describe results in three models: the 2D Ising ferromagnetic model, the 3D Vicsek flocking model and a small-world neuronal network model. We illustrate how the feedback of the autocorrelation function of the order parameter fluctuations is able to shift the system towards its critical point. Since the results rely on universal properties they are expected to be relevant to a variety of other settings.

Equations

$$C(k,t)=C_0(k)g\left(\frac{t}{{\tau }_0(k,\xi )};k\xi \right),$$

$${\tau }_0(k,\xi )=k^{-z}\Gamma (k\xi ;\xi /L)={\xi }^z\Omega (k\xi ),$$

$$\frac{C(k,t=\delta t)}{C(k,0)}$$

$$\frac{C(k=0,\delta t)}{C(k=0,t=0)}\sim 1-A{\xi }^{-z}\sim 1-A(T-T_c{)}^{z\nu },$$

$$d_i=d_{i-1}\text{sign}[AC(1{)}_{(i)}-AC(1{)}_{(i-1)}],$$

$${\delta }_i=(1-AC(1{)}_{(i)}{)}^2,$$

$$T_{(i+1)}=T_{(i)}+\delta \ast d\ast \kappa ,$$

$${\mathbf{v}}_{\mathbf{i}}(t+\Delta t)$$

$${\mathbf{r}}_{\mathbf{i}}(t+\Delta t)$$

$$\phi \equiv \frac{1}{N{v}_0}\left|\sum _{i=1}^N{\mathbf{v}}_i\right|.$$

Full text

The last decade has witnessed an escalating interest in complex biological phenomena at all levels including macroevoluction, neuroscience at different scales, and molecular biology. The observed complexity in nature is often traced to critical phenomena because it resembles the complexity found for critical dynamics in models and theory bak ; mora ; beggs ; reviewBiophys ; chialvo ; tang ; miguel ; flocks . More specifically, it seems that many biological systems reach a “sweet spot” where they attain maximal susceptibility, i.e., sensitivity to changes in the environment, while maintaining internal order.

At present it is not clear how such a critical state can be reached or even maintained. For a complex system like the brain, one might imagine that its control parameters be hard-wired genetically, selected by a long evolutionary process to at a critical point that is biologically most advantageous for survival. However, the critical values of the control parameters depend on system size DOMB , and thus for biological systems to take advantage of critical dynamics they would need to adjust the control parameter as systems contract or expand. We exemplify this problem in Fig. 1 (inset) which sketches how the peak of the susceptibility, the property to be maximized, shifts as the system get larger.

While some physical systems might be large enough that one can assume they are asymptotically near the thermodynamic limit, we note that most biological systems are of moderate size, and finite-size effects are in principle to be expected cavagna2014 . Hence, if the critical point is the best (or only) functioning state for a given biological system, in order to attain it, Darwinian evolution instead of furnishing a set of specific values for the control parameter must allow for a control mechanism such that systems can reach and stay close to a critical point. For a control mechanism to be biologically plausible it should utilize only information that either is completely global or completely local. Here we explore a possibility for such a mechanism.

We show that the first autocorrelation coefficient $AC(1)$ of the order parameter fluctuations can be used to tune a system to the vicinity of its critical point. This is possible because $AC(1)$ peaks at the same point as the susceptibility, yet does so more smoothly than the susceptibility.

This can be understood from dynamic scaling. The dynamic scaling form of the time correlation is HHREVIEW

[TABLE]

where $k$ is the observation wavevector, $\xi$ is the correlation length, the function $g$ is such that $g(t=0)=1$, i.e. $C_{0}(k)$ is the static correlation function, and the characteristic time obeys

[TABLE]

where $g$, $\Gamma$ and $\Omega$ are unspecified scaling functions and $z$ is the dynamic scaling exponent. Now

[TABLE]

For a global quantity, $k=0$ (See Suplementary Material) so that

[TABLE]

where $A>0$ is a time dependent constant. Hence, the normalized time correlation has a maximum at $T=T_{c}$, for fixed $\delta t$.

The main idea is demonstrated here by applying it to three well understood systems, namely the ferromagnetic Ising model, the Vicsek model of flocking and a typical neuronal small-world network. We remark that the results are general enough to be also expected in many other systems.

Ising model.

Fig. 1 illustrates the typical behaviour of the 2D ferromagnetic Ising model at increasing temperatures. The system undergoes a second order phase transition at a critical temperature $T_{c}$, reflected in a steep change in magnetization as well as a sharp peak in susceptibility (Fig. 1A). Equally distinct changes are also demonstrated for the correlation properties of the model computed from appropriate system variables (Fig. 1B). A sharp increase in the average pairwise correlations is observed as the system approaches $T_{c}$, where the correlation length matches the size of the system. The relatively sharp changes in the spatial correlations contrast with the relatively smoother changes in the temporal correlations, as reflected by the first auto-correlation coefficient $AC(1)$ of the magnetization fluctuations around the mean, which at $T_{c}$ approaches unity.

Now we asses how to control the Ising model to stay at the vicinity of the susceptibility peak. According with the discussion in the introduction, we must restrict ourselves to do it using only either local or global information. In that sense, the time correlations evaluated by $AC(1)$ meet such conditions, because it can be computed from a temporally delayed version of a global average of magnetization. In turn, magnetization can be assessed simply by averaging samples of a relatively large number of sites.

To demonstrate control we proceed by choosing an initial random temperature and simulate the dynamics for some large number of Montecarlo (MC) steps, which we denote as an “adaptive iteration step” indexed by $i$. We proceed by estimating the $AC(1)$ of the fluctuations around the mean magnetization during the lapse of time corresponding to the adaptive iteration step $i$ and monitor the change of $AC(1)$ between two consecutive steps $i$, defining

[TABLE]

so that $d$ changes sign when a decrease in $AC(1)$ is detected. We then use the gradient to its maximum value

[TABLE]

to change the future temperature $T(i+1)$ according to

[TABLE]

where $\kappa$ is a small constant that determines how slowly the temperature is adjusted. Its exact values is not crucial for the present results. Successive iterations of Eqs. 5–7 demonstrate convergence of the temperature to the expected value at equilibrium $T_{c}\sim 2.3$. Fig. 2 illustrates typical results for various initial temperatures, which in all cases converge to the vicinity of $T_{c}$. We note that the successive values of the parameters (order, control and $AC(1)$) obtained during the adaptive simulations over-imposes well (i.e., matches) those obtained from equilibrium simulations (i.e. the data of Fig. 1, see Suppl. Material).

Vicsek model.

We were also able to use the $AC(1)$ function to control the Vicsek model VM , the archetypal model for flocking behavior, towards its critical point. In this model, $N$ self-propelled particles endowed with a fixed speed $v_{0}$ move in $d$-dimensional space. At each time step, positions $\mathbf{r}_{i}(t)$ and velocities $\mathbf{v}_{i}(t)$ are updated according to

[TABLE]

where $S_{i}$ is a sphere of radius $r_{c}$ centered at $\mathbf{r}_{i}(t)$. The operator $R_{\eta}$ normalizes its argument and rotates it randomly within a spherical cone centered at it and spanning a solid angle $\eta\Omega_{d}$, where $\Omega_{d}$ is the area of the unit sphere in $d$ dimensions ($\Omega_{2}=2\pi$, $\Omega_{3}=4\pi$).

The order parameter, which measures the degree of flocking, is the normalized modulus of the average velocity VM ; reviewVicsek ,

[TABLE]

$\varphi\in[0,1]$, with $\varphi=O(1/\sqrt{N})\sim 0$ in the disordered phase and $\varphi=O(1)$ in the the ordered phase. We choose $\Delta t=r_{c}=1$, so that the control parameters are the noise amplitude $\eta$, the speed $v_{0}$ and the number density $\rho=N/V$, where $V=L^{d}$ is the volume of the (periodic) box. We apply Eqs.5–7 to this model, using $\eta$ as control parameter and keeping the density fixed. For comparison we over-plotted results from equilibrium runs with values taken during adaptive control of the simulations (Fig. 3). The close match demonstrates that the technique is able to control the flock model near its critical noise amplitude, $\eta_{c}\sim 0.48$ (see further details in Suppl. Material).

Neuronal network model.

Successful control was further demonstrated for a neural network model haimovici consisting of a network of interconnected nodes together with a dynamical rule. The model exhibits a second order phase transition Mahdi on a region of parameters. The model matrix of interactions follows a small-world topology and each node exhibits discrete state excitable dynamics, following the Greenberg-Hastings model Greenberg . Briefly, each node is assigned one of three states: quiescent $Q$, excited $E$, or refractory $R$, and the transition rules are: 1) $Q\rightarrow E$ with a small probability $r_{1}$ ($\sim 10^{-3}$), or if the sum of the connection weights $w_{ij}$ with the active neighbors ($j$) is higher than a threshold $T_{h}$, i.e., $\sum w_{ij}>T_{h}$ and $Q\rightarrow Q$ otherwise; 2) $E\rightarrow R$ always; 3) $R\rightarrow Q$ with a small probability $r_{2}$ ($\sim 10^{-1}$) delaying the transition from the $R$ to the $Q$ state for some time steps. Parameters $r_{1}$ and $r_{2}$, which determine the time scales of self-excitation and of recovery from the excited state, respectively, were kept fixed and $T_{h}$ was updated according to control Eqs. (5–7). The density of active nodes, i.e. in state E, in each time step was taken as the order parameter. As shown for the previous models, $AC(1)$ was able to move and maintain the system near its critical point (here $T_{h}\sim 0.16$) (Fig. 4).

The present results applies, with some differences, to $1^{st}$ or $2^{nd}$ order phase transitions and also for low dimensional dynamical systems exhibiting continuous or discontinuous bifurcations from fixed points to limit cycles which can be controlled near the bifurcation point (see Suppl. Material).

In conclusion we have demonstrated, in three paradigmatic cases, how the autocorrelation function of the order parameter fluctuations allows to establish a feedback loop for the control parameter to shift the system towards its critical point. Our results build on two previous lines of work which come close to describe the control strategy. One is the view of self-organized criticality bak as a feedback between order and control parameters model1 . The other line relates to forecasting of an upcoming tipping point via the generic slowing down of criticality tipping ; tipping2 . The current results go beyond these previous approaches by demonstrating a mechanism that may explain the presence of criticality in some systems, and furthermore providing a strategy of control amenable of practical implementations in different areas. For instance, in neuroscience, this approach could be realized with optogenetical targeting hollopaper to clamp cortical networks to any desired dynamical state, helping to predict its influence on perceptual performance.

Acknowledgements.

Supported by grant 1U19NS107464-01 from the NIH BRAIN Initiative (USA) and by CONICET (Argentina).