Intermittent dynamics in complex systems driven to depletion

When complex systems are driven to depletion by some external factor, their non-stationary dynamics can present an intermittent behaviour between relative tranquility and burst of activity whose consequences are often catastrophic. To understand and ultimately be able to predict such dynamics, we propose an underlying mechanism based on sharp thresholds of a local generalized energy density that naturally leads to negative feedback. We find a transition from a continuous regime to an intermittent one, in which avalanches can be predicted despite the stochastic nature of the process. This model may have applications in many natural and social complex systems where a rapid depletion of resources or generalized energy drives the dynamics. In particular, we show how this model accurately describes the time evolution and avalanches present in a real social system.


Introduction
Avalanches or bursts of activity are present in many natural and social systems that range from sand piles 1 , earthquakes 2 , solar flares 3 , plastic deformation 4 , domain flips in ferromagnets 5 , species mass extinctions 6 , neuronal activity 7 , forest fires 8 , landslides 9 , clustering of self driven particles 10 and hydrodynamics 11 , to daily website hits 12 , attendance to motion pictures 13 , financial markets 14 and book sales 15 .Understanding and ultimately predicting the dynamics of these avalanches is of great scientific and practical importance, as they can lead to catastrophic events, of which large earthquakes and stock market crashes serve as two dramatic examples.
A very fruitful approach to understand the dynamics of many of these systems has been to assume local and relatively simple interactions of their constituting elements which drive the system to a state that displays emergent behaviour.Importantly, this emergence results without the need of an external tuning parameter 16,17 .This is the premise of Self Organized Criticality (SOC) 1 , for which models now abound.A paradigmatic example is that of earthquakes, which may be modelled as dissipative systems presenting SOC 18,19 .There, a slow energy input eventually leads to relatively fast breakdown events spanning many orders of magnitude.Variations of this model 20 , based on the Bak, Tang and Wiesenfeld (BTW) sand-pile model 21 , are able to recover the main characteristics of earthquakes statistics.Roughly speaking, this class of BTW models have two main ingredients.First, there is a continuous input of energy into the system that may continue indefinitely.Second, bursts of activity occur due to excess stress accumulation which is then toppled to neighbouring constituents giving rise to avalanches. Figure 1: Schematic diagram of the negative feedback mechanism proposed: A decaying input signal q(t) is randomly distributed at time t among N (t) elements.When T i (t) crosses a threshold x i , element i becomes inactive, thus decreasing the overall number of active elements N (t).To exemplify this mechanism, in the upper right panel elements labeled i = 2, 4, 5, 8 and 9 become inactive at time t + 1 since the condition T i (t) < x i was fulfilled for each of these i s at time t.
There exist, however, a plethora of other systems which are instead characterized by a steady decrease of their activity all the way to complete inertness.This happens as a result of a steady decay of the input energy, which translates into the system as an effective stress to which it must adapt, reaching at every step a new equilibrium configuration.While in some instances this is achieved in a relatively smooth way, in others the system is unable to adapt fast enough giving rise to large, often catastrophic events.Indeed, one may see an intermittent behaviour between relative tranquility and bursts of activity.Actually, it is now recognized that continuous driving forces can induce transitions between continuous and intermittent dynamics 22 .This feature is of great importance, since systems that seem to be adapting well (smoothly) to such external stresses, may suddenly break down without any apparent warning signs.Examples of systems presenting this general behaviour range from the group of theaters showing a particular movie as weekly audience decreases, the number of surviving a mammalian genera during an extinction period or the price of a stock during a sell-off.

Model definitions
To model this behaviour we consider a system of non-interacting elements that receive (but do not accumulate) a share of the energy externally provided.Each element is pre-assigned a sharp threshold such that if the energy it has received is below this threshold, the element is permanently removed from the system.These thresholds are considered as intrinsic properties of the elements, and therefore do not change with time.We show below that this simple model naturally leads to a negative feedback responsible for the adaptation of the system as a whole.The general mechanism proposed is depicted in figure 1.This model reproduces both the smooth response as well as the intermittent one and provides a means of predicting catastrophic events in complex systems as they are driven to extinction.To be mathematically precise, our model consists of N basic elements labelled with i = 1, . . ., N .We then introduce a time-decaying external input signal q(t), that, for sake of clarity, we assume takes integer values and decays rapidly with time.This function represents the total generalized available energy, which is randomly distributed among the N (t) active elements at each time step t.An element belonging to this group becomes -and remainsinactive when the amount of energy allotted to it is below its own, pre-assigned threshold.To achieve this, a fixed set of time-independent thresholds x i > 0 for each of the constituent elements i = 1, . . ., N is drawn randomly from a distribution.We stress the point that at any given time step each element can be either active or inactive, but once it has become inactive it remains in that state.Let σ i (t) be the state of element i, with σ i (t) = 0 for the element being inactive, or one otherwise.If N (t) = N i=1 σ i (t) denotes the number of active elements at time t, the system evolves according to the following simple microscopic dynamical rule: Here Θ(x) represents the Heaviside step function, and T i (t) is a portion of the input signal uniformly and randomly allotted only to active element i such that i=1 T i (t) = q(t).Henceforth, we will call the set of variables T i (t) the occupancies.Notice that at a given time t, the probability of finding a particular configuration {T 1 (t), . . ., T N (t) (t)} follows a multinomial distribution with uniform probability 1/N (t).This implies that the average occupancy is the same for all elements and equal to We take the initial conditions in which all constituents are active at t = 0, that is N (t = 0) = N , and define the activity for this model as the number of elements removed at each step ∆N (t) = One can gain some intuition on how the system behaves under this mechanism by looking at the time evolution of the density of active thresholds and the density of occupancies, denoted as ρ t (x) and ρ t (T ) respectively.This will become useful to gain a good grasp on the mean-field solution to the model, presented below.Figure 2a shows a sketch of the dynamical evolution of the system in terms of these two densities.Let us assume that at a given time these two distributions do not overlap (Figure 2a, t = 0).Then, as time evolves, the distribution ρ t (T ) moves to its left since the average occupancy T (t) = q(t)/N (t) decreases over time, as q(t) is a decreasing function and N (t) has not changed yet.Eventually, the two distributions do overlap (Figure 2a, t = 1).
This means that the overlapping fraction of active thresholds will become inactive, increasing the average occupancy T (t), and pushing, in turn, the distribution ρ t (T ) back towards higher values.
As a consequence, a non-zero activity indirectly hinders any further activity, in what constitutes a negative feedback mechanism.Now again, the distributions are not overlapping (Figure 2a, t = 2) but, as time evolves, this mechanism reoccurs (Figure 2a at t = 3 and t = 4).This sequence of events repeats itself until all elements have become inactive.Since in this sequence the system goes through periods of zero activity, Figure 2a actually illustrates its typical response when it is in the herein called intermittent regime.But with this representation it is also possible to foresee that if the external signal decays relatively slowly, the system will be able to adapt to the depletion of resources, so that small but finite portions of ρ t (x) are removed at every iteration (Figure 2c).We call this the continuous regime since the activity is always non-zero.The different behaviour of the system in both regimes can be appreciated in the plots of N (t) and ∆N (t) presented on Figures 2b and 2d for the intermittent and continuous regimes, respectively.These plots also help illustrate an important characteristic of the model: the maximum activity can be larger in the intermittent regime by over an order of magnitude as compared to the continuous one.
At this point it is important to contrast our model with the BTW model, as both present avalanches and are based on a sharp threshold dynamics as a response to some external energy input.In the BTW model, the elements q i (t) may represent the accumulated stress allocated in one of the N (t) geological faults, and, as in our model, an earthquake or avalanche occurs when a fault crosses its threshold.But, unlike our model in which the elements remain inactive once they are brought to that state, as soon as a geological fault has toppled its excess stress, it can start again accumulating new stress.Another important difference is that in the BTW model the input signal keeps providing energy which the system continuously dissipates, while, in our case, the input signals evolves in times towards zero.Furthermore, our model does not allow any stress accumulation; in other words, it is history independent.Indeed, despite their shared characteristics, the actual mechanism giving raise to bursts of activity in each model is completely opposite in nature: in the BTW model avalanches result from a positive feedback mechanism, while in ours, a negative feedback one plays the central role.

Results
For concreteness, to quantify the behaviour of the model we have chosen an exponentially decaying external signal q(t) = q 0 e −t/τ and proceeded to identify the continuous and intermittent regimes of the model in its space of parameters.As relevant parameters, apart of τ , we have also chosen the initial coefficient of variation κ ≡ σ/µ of the distribution of active thresholds ρ t=0 (x), where µ and σ 2 are its first two cumulants.We have subsequently estimated the dynamical phase diagram corresponding to the the microscopic dynamics of equation (1) using two methods: i) by deriving a simple dynamical mean-field equation and ii) by Monte Carlo simulations.
Even though it is possible to obtain an exact effective macroscopic dynamics using the method of the generating functional analysis a la Martin-Siggia-Rose 23 , the resulting equations are rather cumbersome to analyse.Thus we have opted to derive a simpler dynamical mean-field equation, which turns out to describe the exact microscopic dynamics remarkably well.After some algebra, a dynamical mean-field equation for the average number of active elements, which we also denote as N (t), is given by (see SI) To identify in which regime the model is in the parameter space (κ, 1/τ ), we introduce two order parameters.The first one is a natural way to quantify intermittency through the normalized cumulative time of inactivity (see SI).In particular, this quantity is 0 when the system is in the continuous regime, and non-zero when in the intermittent one.The second parameter characterizes the ability of the system to follow the input signal and corresponds to see whether z(t) equals a constant z after some transient.By iterating the dynamical mean field equation ( 2) and following the values of these two parameters we are able to identify four regimes or phases in the (κ, 1/τ )-plane.Moreover, we are able to obtain analytical expressions for some of the transition lines.These four regimes are shown in Figure 3a.The intermittent regime, depicted in Figure 3a by a red filled area, is characterised by the impossibility of the system to synchronise with the external signal, resulting in avalanches.The continuous & asynchronous regime I (shown as the dark blue area in the same figure) corresponds to the system being again unable to synchronize to the external signal, but this time in a smooth way, i.e., avalanches are no longer generated.On the other side, in the continuous & synchronous regime, shown by the green area in Figure 3a, the system is able to follow the external signal.Finally, the continuous & synchronous regime II, captures the impossibility of the system to match the external signal due solely to restrictions in the initial conditions T (t) > 0.
Notice that, within the mean-field approximation, it is possible to obtain analytically the bifurcation line, denoted as κ c (τ ), separating the asynchronous to the synchronous regimes, as well as the bifurcation line κ na (τ ) separating the synchronous to the asynchronous regime II.The first line is derived by assuming that the synchronous regime is characterised by z(t) = z being a constant.This implies, in turn, that z and τ are related by the formula z (τ ) = erfc −1 (2e −1/τ ).
Naturally, the transition to the intermittent regime is derived by assuming a continuous bifurcation of z(t) around z .After some algebra, the line separating both regimes in the (κ, 1/τ )-plane The Monte Carlo simulations were carried out by taking q(0) = 2 × 10 9 , N (0) = 2 × 10 5 , and µ = 1000, and subsequently varying the parameters τ and σ, to identify the regimes in the phase diagram.Results of 50 different initial random distributions for ρ t (x) were averaged.For the mean field model, the system was considered to be inactive whenever (N (t)−N (t+1))/N (t) < 1×10 −5 .
is given by the following set of coupled equations: with 1 τ = − log erfc(z )

2
. Thus, given µ, the solution of (3) results in a pair (κ c (τ ), T (τ )), where the first function κ c (τ ) corresponds to the bifurcation line separating the aforementioned regimes, while line T (τ ) tells which value T (t) takes precisely at the transition.As shown in Figure 3a, the agreement between the bifurcation line (3) and the dynamical phase diagram obtained from the same equation is excellent (see the SI for more details on the synchronization plots for this and other values of µ).The second bifurcation line, as reported in Figure 3a, follows from the fact that in the continuous and synchronous regime we must have and, since N (t) decays precisely like the signal q(t), one can easily derive from equation ( 2) that . However, as T ≥ 0 there is a part of the phase diagram inaccessible in the synchronous regime, given by the following line κ na (τ ) in the (κ, 1/τ )-plane (see SI): How well the mean field description captures the microscopic model can be appreciated by comparing Figures 3b and 3c, where we have calculated the cumulative time of inactivity for each model, using Monte Carlo simulations for the microscopic one.Despite the simplicity of equation ( 2), the agreement with the microscopic model is reasonably good.Notice that these plots show there is a clear border separating the intermittent and synchronous regimes, and that, in general, intermittency can be avoided for faster decays with relatively wide ρ t (x) distributions, as expected.
Furthermore, both models present a gap: below a critical value of 1/τ , the system is continuous regardless of how small κ is.The magnitude of this gap is actually a function of µ, and its analytical expression can be found for the mean-field description (see SI).We point out that, unlike the meanfield model, the microscopic one does not present a synchronized regime.Rather, in the latter the activity always lags behind the external signal, i.e., it is always marginally desynchronized.We find, however, that the activity does become less synchronized in the intermittent regime, but no clear transition is observed for this model.

Application to a real system
Despite its apparent simplicity, the model presented here is robust enough to be applicable to real complex systems.As a proof of principle, we now implement it to investigate the dynamics of the supply of theaters for a given movie that played in the US from 1970 to 2000.The dynamics of the attendance to movie theatres in the U.S. has been studied recently in Ref. 13 , where it was shown that the attendance is in general very well described by a decaying exponential (see example of left plot in Figure 4) as inferred from the weekly revenue of a movie (see SI), and is thus ideally suited to test our model.This system meets our model if we consider that as q(t) people fill up the N (t) theaters at time t, a theater removes the movie at time t + 1 if the local attendance T i (t) is lower than its threshold x i .In other words, a theater removes a movie after the revenue from one week is below some pre-established value.We assume now that the distribution of local thresholds is again normally distributed.As a first order approximation we assume that at time t the whole available population q(t) is randomly distributed over the available theaters N .We realize that was obtained from www.boxofficemojo.com.See Ref. 13 for details.
this is an enormous simplification of the actual problem since, of course, geographic restrictions limit the available theaters for people.Nevertheless, similar results as the ones presented below are obtained by sectioning the data and assigning smaller variance per section.
The right panel in Figure 4 presents four examples of how our model (thick pink line and square symbols) is capable of reproducing the theater dynamics (black dots) remarkably well, in some cases spanning almost three orders of magnitude.Note the logarithmic scale on the number of theater axis of this figure.Both small and relatively large events are equally reproducible.In all examples shown, the same parameters µ = 400, σ = 0.45µ, were used.In other words, we only require the initial number of theaters to reproduce the whole time series of N (t) using as input the specific time series for q(t) and a normal distribution of thresholds.Other than the fixed choice of distribution of thresholds, there are no free parameters in this model.
We now apply the model to a collection of 3000 different movies that played in the U.S from 1970 to 2000 (see SI).We postulate that this social system is characterized by a single distribution of thresholds, and therefore, the same parameters µ = 400, σ = 0.45µ are used in all instances, and only q(t) and N (0) are different for each movie and determined by external factors (see SI for further details).To assess the capability of the model to completely recover the dynamics of this social system we obtain the distribution of ∆N (t) normalized by N (0) for each movie from both the model (solid line in figure 5) and the actual data (dots in figure 5).Notice that the model recovers the whole shape of the distribution of these normalized avalanches.This distribution is actually well fit to a power law with an exponential tail (not shown).In this particular system, an epidemic branching mechanism is responsible for the exponential decay of the effective energy q(t), Ref. 13 .It is possible that other complex systems also present a similar rapid decline of resources, driving the dynamics towards intermittency.

Conclusions
To conclude, we have introduced a negative feedback model with sharp thresholds to describe the dynamics of complex systems driven to depletion, which is simple enough to be amenable to analytical treatment and yet powerful enough to accurately predict the dynamics in real complex systems.We foresee that this model may be used to explain and predict small and large avalanches in other systems for which there is an exponentially fast decline of some quantity playing the role of an energy source randomly distributed among elements that need a certain minimum amount of it to remain active (or survive).Since large avalanches are indeed catastrophic in many social and natural settings, implementing variations of the model presented in this work could help prevent them.
Supplementary Information is linked to the online version of the paper at www.nature.com/natureAvalanche magnitudes were normalized by the corresponding N (0).A total of 100 runs per movie were averaged to cancel out any special condition of the initial threshold distributions.The distributions ρ t (x) were chosen randomly from a normal distribution with fixed (µ = 400, σ = 0.45) for all movies.Note that the model reproduces both small and large events.

Figure 2 :
Figure 2: Schematics of the time evolution of the number of active elements N (t), represented

Figure 3 :
Figure 3: Analysis of dynamical phase-diagram.a) Phase diagram for the mean field model.Plots

Figure 4 :
Figure 4: Application of the model to real data.Left plot: Attendance to the movie The Dark

Figure 5 :
Figure 5: Avalanches in this complex system.Probability density of the avalanche magnitude for