Time-series thresholding and the definition of avalanche size

TL;DR

This paper reveals that thresholding in continuous time series can obscure true avalanche scaling behavior, leading to potential misinterpretations of neural activity data and false universality class assignments.

Contribution

It demonstrates analytically and computationally how threshold choices affect avalanche size measurements and scaling analysis in continuous stochastic processes.

Findings

Large thresholds can mask true power-law scaling.

Transient effects dominate when using area-based avalanche sizes.

Misinterpretations can lead to incorrect universality class identification.

Abstract

Avalanches whose sizes and durations are distributed as power laws appear in many contexts. Here, we show that there is a hidden peril in thresholding continuous times series --either from empirical or synthetic data-- for the detection of avalanches. In particular, we consider two possible alternative definitions of avalanche size used e.g. in the empirical determination of avalanche exponents in the analysis of neural-activity data. By performing analytical and computational studies of an Ornstein-Uhlenbeck process (taken as a guiding example) we show that if (i) relatively large threshold values are employed to determine the beginning and ending of avalanches and if (ii) --as sometimes done in the literature-- avalanche sizes are defined as the total area (above zero) of the avalanche, then true asymptotic scaling behavior is not seen, instead the observations are dominated by transient effects. These can induce misinterpretations of the resulting scaling regimes as well as to a wrong assignation of universality classes.

Tables

Table 1. Table 1: Summary of the avalanche (mean-field) exponents: size ( τ ) 𝜏 (\tau) , duration ( α ) 𝛼 (\alpha) and averaged avalanche size ( γ ) 𝛾 (\gamma) for the (un-biased) branching process (BP) and the (un-biased) random walk (RW); see e.g. di Santo  et al. ( 2017 ) .

	$P\left(S\right)\sim S^{-\tau }$	$P\left(T\right)\sim S^{-\alpha }$	$P\left(S\mid T\right)\sim T^{\gamma }$
BP	$\tau =3/2$	$\alpha =2$	$\gamma =2$
RW	$\tau =4/3$	$\alpha =3/2$	$\gamma =3/2$

Equations

$$\dot{x}(t)=\sigma \eta (t),$$

$$\dot{x}(t)=-ax(t)+\sigma \eta (t),$$

$$\Sigma (T)=S(T)+s^{\ast }(T)=\tilde{c}{T}^{3/2}+\theta T,$$

$$d\Sigma =(c{T}^{1/2}+\theta )dT$$

$$P(\Sigma (T))=P(T)\frac{dT}{d\Sigma }=\mathcal{N}\frac{T^{-3/2}}{c{T}^{1/2}+\theta }.$$

$$P(\Sigma )\approx \frac{\mathcal{N}}{c}{T}^{-2}\approx \frac{\mathcal{N}}{c}[(\Sigma /c{)}^{2/3}{]}^{-2}\sim {\Sigma }^{-4/3},$$

$$P(\Sigma )\approx \frac{\mathcal{N}}{\theta }{T}^{-3/2}\approx \frac{\mathcal{N}}{\theta }[\Sigma /\theta {]}^{-3/2}\sim {\Sigma }^{-3/2},$$

Full text

I Introduction

Episodic outbursts of activity or “avalanches” of highly variable durations and sizes are observed in a large variety of scenarios in condensed matter physics (vortices of type II superconductors Altshuler and Johansen (2004) and Barkhaussen noise Alessandro et al. (1990); Cizeau et al. (1997)), high-energy astrophysics (X-ray flares Wang and Dai (2013)), geophysics (earthquakes Bak et al. (2002)), meteorology (rainfall Peters et al. (2001)), neuroscience (neuronal avalanches Beggs and Plenz (2003)), as well as in other biological systems (gene knock-out cascades) Roli et al. (2015)) and man-made systems (failures on electrical power grids Kinney et al. (2005)). The probability distributions of sizes and durations of such avalanches often exhibit a “fat-tail” that can be fitted as a power-law distribution; i.e. the fingerprint of scaling behavior. Such scaling or scale invariance is often considered as evidence of underlying criticality and many of the above systems are claimed to operate at (tuned or self-organized) critical points Bak (2013); Bak et al. (1987); Jensen (1998); Rubin et al. (2014); Dickman et al. (2000). In particular, in the context of biology the idea that living systems (parts, aspects or groups of them) may extract important functional advantages from operating at criticality –i.e. at the edge of two different phases– has been deeply explored in recent years Mora and Bialek (2011); Muñoz (2018).

In this regard, groundbreaking experimental evidence by Beggs and Plenz Beggs and Plenz (2003), revealed the existence of scale-invariant episodes of electrochemical activity in neural tissues thereafter named neural avalanches. Subsequently, neural avalanches were robustly detected across a large variety of experimental settings, tissues and species Beggs and Plenz (2003); Plenz and Thiagarajan (2007); Beggs (2008); Lombardi et al. (2012); Petermann et al. (2009); Hahn et al. (2010); Palva et al. (2013); Bellay et al. (2015). In particular, neuronal avalanche sizes, $S$, were robustly observed to be distributed as a power-law $P(S)\sim S^{-\tau}$ with $\tau\approx 3/2$ up to some upper cut-off; similarly, avalanche durations $T$ were well fitted by $P(T)\sim T^{-\alpha}$ with $\alpha\approx 2$ up to some characteristic maximum time Beggs and Plenz (2003). Furthermore, fundamental scaling relationships Muñoz et al. (1999) were observed to be fulfilled: e.g. the averaged avalanche size scales as $\langle S\rangle\sim T^{\gamma}$ and the set of exponents obey $\gamma=(\alpha-1)/(\tau-1)$ Friedman et al. (2012).

This set of empirically reported exponent values is in agreement with that of the well-known critical (or “unbiased”) branching processes, also called Galton-Watson process, originally introduced to describe the statistics of the extinction of family names) Watson and Galton (1875); Liggett (2004); Zapperi et al. (1995); Redner (2001). Actually, the set of exponent values $\tau=3/2$, $\alpha=2$ and $\gamma=2$ are extremely universal as they are shared by many different propagation processes in high-dimensional systems as well as in many types of networks Alessandro et al. (1990); Karrer et al. (2014). In particular, they are the mean-field exponents shared by models such as the contact process, directed and isotropic percolation, susceptible-infected-susceptible, and a large list of other prototypical models for spreading/propagation dynamics above their respective upper critical dimensions Muñoz et al. (1999); Binney et al. (1993); Dobrinevski et al. (2015); Garcia-Millan et al. (2015); Bonachela et al. (2010).

Thus, it was conjectured that neuronal systems might operate close to the edge of marginal propagation of neural (electro-chemical) activity Beggs and Plenz (2003); Chialvo (2004), opening the door to exciting theoretical perspectives and some debate (see Muñoz (2018) for a recent review). However, as extensively discussed in the literature, diverse generative processes for the emergence of power-laws exist Beggs and Timme (2012); Newman (2005); Mitzenmacher (2004), and not all power-law distributions can be taken as a signature of criticality. For instance, a diverging correlation length needs to be identified to assign a given phenomenon to criticality. In recent years, some authors have suggested that the origin of the observed power-law scaling in neural systems might stem from other types of criticality (rather than marginal propagation) Martinello et al. (2017); di Santo et al. (2018, 2016) or even be unrelated to critical behavior Touboul and Destexhe (2017); Priesemann and Shriki (2018).

In the present brief paper –leaving aside the putative connection with criticality– we contribute with an additional piece of information to the already controversial discussion about the statistics of neuronal avalanches. In particular, we show that some of the reported empirical evidence in favor of the value $\tau=3/2$ –and thus, seemingly in favor of the existence of an underlying critical branching process– might be misleading as there is a technical problem in the way avalanches are measured, which hinders the observation of the true asymptotic behavior. More in general, we underline that particular attention needs to be taken when avalanches of activity –defined by thresholding– are inferred from a continuous time series of activity. Our findings, add to the recent literature warning on the “perils” associated with thresholding in timeseries 111For instance, since the same avalanche could e.g. be split in two by the effect of a higher threshold, it would introduce correlations between different avalanches, altering the measurement of exponents Laurson et al. (2009); Font-Clos et al. (2015). Laurson et al. (2009); Font-Clos et al. (2015).

II Definition of avalanches

As discussed in the introduction, avalanching phenomena are best described, at least in mean-field, as branching processes. Such processes can reach the value [math] which is an absorbing state: avalanches are naturally defined as excursions away from such a state caused by small perturbations di Santo et al. (2017). However, in many contexts –including neuroscience but not only– the term “avalanche” is used to refer to excursions of time series above some given (arbitrary) threshold, regardless of absorbing states the existence of any absorbing state. In this section we discuss such avalanches and their statistics.

II.0.1 Avalanches in the Wiener and Ornstein-Uhlenbeck processes

Let us consider, for argument’s sake, a time series for a stochastic real variable $x$ –as illustrated in Fig.1– generated by a Wiener Process, i.e. by a continuous-time unbiased random walk (RW) defined by the following Langevin equation Gardiner (2004):

[TABLE]

where $\eta(t)$ is a Gaussian white noise with zero mean and unit variance, and $\sigma$ is the noise amplitude. For such a time series (which can be thought as describing the time course of the activity of some arbitrary system) the duration $T$ of an avalanche is the amount of time for which $x$ stays above a given threshold, i.e. an avalanche begins/ends when the activity signal crosses the threshold from below/above; the avalanche size $S$ is the area covered between the walk trajectory and the threshold reference line (see Fig. 1). Observe that similarly, given the symmetry of the process, one could also define avalanches as excursions below threshold.

The probability distribution of avalanche durations $T$ can be straightforwardly identified with the first-return time statistics of random walks (see Fig.1) which is well-known to scale with an exponent $\alpha=3/2$. Similarly, also the size-distribution exponent $\tau=4/3$ and the remaining exponent $\gamma=3/2$ are well-known for random walks (see Fig. 1 and Table I); pedagogical derivations of these results, as well as a comparison with the branching process class can be found in e.g. di Santo et al. (2017); Artime et al. (2018). Importantly, these results for the random walk do not depend of the value of the chosen threshold.

In the more general case in which the walker is confined to hover around a given mean value, one can describe the problem, in first approximation, as an Ornstein-Uhlenbeck process Gardiner (2004):

[TABLE]

where there is an additional linear force term, $-ax$ (corresponding to the negative derivative of the parabolic potential bounding the walker close to $x=0$). Such a force introduces an upper cutoff in the first-return times statistics of the unbiased RW (see e.g. Artime et al. (2018) for a detailed derivation). Thus, avalanches intended as excursions above a given threshold in a process with a well-defined steady-state value, have power-law-distributed sizes and durations, with the exponents of the RW class (as in table I and as in Fig.1), but only up to an upper-cut-off scale controlled by $1/a$, such that it goes to infinity, i.e. it disappears when $a$ vanishes.

As a corollary of all this, let us remark that many real time series describing (e.g. biological) problems in which some stochastic variable fluctuates symmetrically around a given mean value exhibit effective avalanching behavior that –up to certain scale of size and time– can be described by the exponent values of the RW. Let us stress again that, as argued above, it can be a matter of debated whether this type of behavior –describable in terms of random-walk excursions above a threshold– can be called “avalanching”. Actually, for most of the examples of interest in physics, as discussed in the first paragraph of the Introduction, this does not constitute an adequate description as it does not include any absorbing state.

II.0.2 On the definition of avalanche size

In many other circumstances time series exhibit *asymmetric * excursions around a mean value and/or may become trapped at some absorbing state. An example of this is obtained when the variable under scrutiny is a positive-definite density (e.g. of neural activity), which, by definition, is constrained to take positive values, $x(t)>0$. In such cases, especially in the ones when there is always some lingering activity so that the zero-value is hardly reached (see Fig.2) a threshold $\theta>0$ is often employed to define avalanches as periods during which the activity remains above such a threshold. In these situations, two alternative possibilities are often used in the literature to measure the size of a so-defined avalanches Dalla Porta and Copelli (2018):

(A) Following the random walk analogy, as done above, for a given avalanche, one can define its size $S$ as the area in between the time-series curve and the threshold ($\theta$) reference line.

(B) Alternatively, one can define the avalanche size $\Sigma$ as the overall integral of the time series during the avalanche, i.e. above the reference line $x=0$ (see e.g. Poil et al. (2012); Larremore et al. (2014) but there are other works making this choice).

The difference between the two criteria to define avalanche sizes is sketched in Figure 2. $\Sigma$, the total integral of the activity is equal to $\Sigma=s^{*}+S$, where $S$ is the integral of the signal above threshold, and $s^{*}$ is the area of the rectangle under the threshold.

In what follows, we compare the statistics of avalanches obtained using these two alternative definitions of size A and B for an Ornstein-Uhlenbeck process. This will serve as an illustration of a more general phenomenon that may also occur for other processes, such as the one sketched in Fig.2.

First we discuss computational results and then we employ scaling arguments to explain the findings. On the one hand, as already shown in Fig.1 using $S$, i.e. criterion A, one reproduces the expected theoretical results for all three avalanche exponents. On the other hand, as illustrated in Figure 3, the statistics of avalanche sizes, as determined employing $\Sigma$ for an Ornstein-Uhlenbeck process is anomalous and does not match the expectations for the theoretically known values i.e. the measured value $\tau\approx 3/2$ does not coincide with the expected value $\tau=4/3$. In particular, the numerical observation of, $\tau\approx 3/2$ could be (wrongly) taken as evidence of branching process-like scaling Poil et al. (2012); Larremore et al. (2014). The fact that there is something suspicious with the definition (ii) can be noticed observing that both sizes and durations scale in the same way, entailing $\gamma=(\alpha-1)/(\tau-1)=1$, which would imply a locally linear (i.e. “tent like”) shape of avalanches Baldassarri et al. (2003); Papanikolaou et al. (2011); Laurson et al. (2013).

II.0.3 Scaling arguments

The correction $s^{*}$ for a given avalanche (such that $\Sigma=S+s^{*}$) is nothing but $s^{*}=\theta T$, where $T$ is the avalanche duration. The distribution of first-passage times for the Wiener process is given by $P(T)\sim T^{-3/2}$. As said above, the same result holds for the Ornstein-Uhlenbeck case up to an upper cut-off. Thus, $\Sigma$ has a correction $s^{*}$ with respect to $S$ that scales as the avalanche duration: $P(s^{*})\sim{s^{*}}^{-3/2}$. Assuming that the probability to observe a given size $S$ conditioned to a given avalanche duration $T$, $P(S|T)$, is a peaked function around its mean value (as usually occurs for avalanches Muñoz et al. (1999)) and using the fact that $\langle S\rangle\sim T^{\gamma}$, with $\gamma=3/2$, then

[TABLE]

from where it follows that

[TABLE]

Thus, we can readily write (using the implicit function theorem):

[TABLE]

From this, in the limit of vanishing threshold $\theta$ in Eq.(4), one has

[TABLE]

which is the correct result for the avalanche size distribution of an Ornstein-Uhlenbeck process. On the other hand, for larger values of the threshold $\theta$ and relatively small values of $T$ (and, thus, also typically small values of $\Sigma$) one has

[TABLE]

in agreement with the numerical observation above. In other words, the additional contribution $s^{*}$ dominates the scaling behavior of the avalanche size $\Sigma$ when $\theta$ is relatively large. It is important to emphasize that, in any case, one should recover the correct asymptotic value –i.e. the behavior for large values of $\Sigma$– of the avalanche size exponent ($\tau=4/3$) for any value of $\theta$ but this requires going to larger and larger avalanche sizes as $\theta$ is chosen larger and larger. In particular, Figure 4 illustrates that there is a crossover from the value $\tau\approx 3/2$ measured for small avalanche-sizes to the true asymptotic scaling $\tau=4/3$, for larger sizes. The crossover point grows with $\theta$, so that the effect is not observed for $\theta\approx 0$, but may extend for many scales even for moderate values of $\theta$. In particular, given that an upper cut-off to scaling may exist (controlled e.g. by $1/a$ in the case of an Ornstein-Uhlenbeck process that we are considering here or by finite-size effects), the transient behavior usually extends all the way up to the cut-off, so that the true asymptotic behavior can be unobservable with criterion B if large values of $\theta$ are considered.

Thus, summing up, considering criterion B for the definition of avalanche sizes together with relatively large threshold values or not sufficiently large statistics, may lead to the observation of an effective value $\tau\approx 3/2$; this may induce a misinterpretation of the scaling universality class, suggesting it is branching-process-like rather than what it actually is: a random-walk-like process.

Let us emphasize that the previous discussion has been done for an Ornstein-Uhlenbeck process. However, it perfectly illustrates the problem associated with criterion B in more general circumstances, e.g. it also applies to asymmetric processes as the one sketched in Fig.2. In any case, criterion B mixes the scalings of actual sizes and times, leading to potential interpretation errors in general stochastic processes. This problem is at the origin of miss-classification of scaling behavior in existing works analyzing neuronal avalanches (see e.g. Poil et al. (2012); Larremore et al. (2014)).

III Conclusion

In this brief paper we have shown that an inappropriate definition of avalanche sizes as measured as excursions above a given threshold in continuous timeseries can lead to misleading conclusions. To illustrate this, we have studied a simple Ornstein-Uhlenbeck process (representing e.g. the time course of activity in a mesoscopic model of neural activity) and have measured avalanches sizes in two possible ways: (i) as the integrated activity $S$ over a given threshold and (ii) integrating the total activity signal in between two threshold crossings, as illustrated in Fig. 2. We have shown both computationally and using scaling arguments that this latest definition can induce strong biases in the determination of the avalanche-size exponent $\tau$.

In particular, if large values of the threshold $\theta$ are considered, then –for relatively small avalanches– one observes the exponent value $\tau\approx 3/2$ which could lead to the erroneous interpretation that an effective un-biased branching process dynamics exists. On the other hand, for sufficiently small threshold values and for sufficiently large avalanche sizes the correct scaling $\tau=4/3$ is recovered. As discussed above the problem associated with criterion B extends to any type of stochastic process as it mixes up the scaling of actual sizes with that of durations, giving rise to misleading results.

This is the underlying reason why recent analyses of avalanches in mesoscopic models of neural activity that consider relatively large thresholds Poil et al. (2012); Larremore et al. (2014); Dalla Porta and Copelli (2018)) obtain $\tau\approx 3/2$ –compatible with the scaling of a critical branching process–, a result that our analysis reveal that is not asymptotic. Their corresponding underlying dynamics describe fluctuations around a given mean value and, thus, the associated avalanches should be related to excursions of random walkers and not to critical branching processes.

As an important final remark, let us stress that it is essential –and it should always be done– to consider the full set of avalanche exponents i.e. $\tau$, $\alpha$, and $\gamma$ as well as the scaling relations between them, in order to avoid possible errors and misleading interpretations and properly identify the type of scaling behavior.

Acknowledgements.

We acknowledge the Spanish Ministry and Agencia Estatal de investigación (AEI) through grant FIS2017-84256-P (FEDER funds) for financial support. This study has been partially financed by the Consejería de Conocimiento, Investigación y Universidad, Junta de Andalucía and European Regional Development Fund (ERDF), ref. SOMM17/6105/UGR as well as to TEACH IN PARMA for support to M.A. M. We thank V. Buendía for very useful comments and a critical reading of the manuscript.