A combinatorial framework to quantify peak/pit asymmetries in complex dynamics
Abstract
We explore a combinatorial framework which efficiently quantifies the asymmetries between minima and maxima in local fluctuations of time series. We firstly showcase its performance by applying it to a battery of synthetic cases. We find rigorous results on some canonical dynamical models (stochastic processes with and without correlations, chaotic processes) complemented by extensive numerical simulations for a range of processes which indicate that the methodology correctly distinguishes different complex dynamics and outperforms state of the art metrics in several cases. Subsequently, we apply this methodology to realworld problems emerging across several disciplines including cases in neurobiology, finance and climate science. We conclude that differences between the statistics of local maxima and local minima in time series are highly informative of the complex underlying dynamics and a graphtheoretic extraction procedure allows to use these features for statistical learning purposed.
pacs:
I Introduction
A major challenge in studying temporally unfolding natural systems is making sense of data that are noisy, reflect processes at local and global scales, and are likely nonstationary. This is exemplified in diverse domains in the natural sciences where theoreticallydriven research aims to describe how system organization and interaction dynamics relate to produce time series features NTSA (); ecoTSA (). Scientists therefore require timeseries descriptors that are: highly informative while still offering a good compression, relatively robust against instrument induced fluctuations, and (particularly in neurobiological domains) invariant to individual differences that exist between samples.
Particularly important for research across disciplines including neurobiology, finance and earth sciences is that such descriptors differentiate the contribution of oscillatory processes from that of pointlike impulses that are inherently nonoscillatory. This is patently apparent when considering spectral power measures, where both oscillatory processes and regularly spaced pointlike events produce similar signatures. Additionally, the structure of pointlike events could itself reflect different processes that determine local minima and maxima in a time series.
This translates into a need to quantitatively measure the anisotropy of a fluctuating signal in a manner that captures to the difference in the structure of peaks (local maxima) and pits (local minima).
Departing from analytic approaches based on consideration of an entire time series – captured by variants of spectral analysis, entropy, or nonlinear signatures – one that contrasts
structure within local peaks and pits is certainly needed, and we present such an approach here.
To further motivate the aforementioned theoretical challenge, we shall discuss the issue in some concrete disciplines. In neurobiology, and more specifically within research examining the brain’s spontaneous mode of operation during wakeful rest (“resting state processes”), it has long been known that
Bloodoxygenleveldependent (BOLD) time series in the human brain show strong power at low frequencies (Hz) and that such frequencies underlie restingstate
connectivity Cordes01 (). This has produced an entrenched view of natural brain dynamics: these dynamics are taken to reflect slow, synchronous oscillations at low frequencies. However, more recent research is shedding light on the problems with such descriptions. For instance, Petridou13 () identified strong and highly infrequent spontaneously occurring spikelike events within recordings of resting state brain activity. They showed that removing these events from the time series reduced lowfrequency power by as much as in certain brain areas. This removal also reduced correlations within brain networks by as much as . Other work suggests that such infrequent spontaneous events carry much weight in explaining important phenomena related to brain function. In particular, Taglizaucchi et al. Tagliazucchi16 () showed that wellstudied connectivity networks in the human brain can be accurately reconstructed even keeping only of the data in each time series (for related work, see Liu13 ()).
In parallel with the shift towards consideration of the role of rare, relatively extreme events, other work began examining in more detail the specific features of both local maxima and local minima in neurobiological time series. Many neurobiological time series, including those measured by fMRI scanners, do not have a natural minima where the measured signal is zero. From the neurobiological perspective, the need to understand minima and maxima in such cases emerges naturally within any model in which oscillatory patterns or power spectra are not the only source of information on brain activity.
To illustrate, in the auditory system, activity peaks are known to track physical features of the input (e.g., frequency or amplitude), so that external stimuli produce well defined steadystate responses (peaks in activity spectra) that track timevarying features of auditory inputs (e.g., Patel00 ()). However, within the same auditory system, activity pits may be further impacted by other factors, including ”dampening” indcued by visual processing, the level of overall attention paid to the auditory stimuli, or the degree of engagement in memory maintenance. This is a simple example of how local maxima and local minima may provide information about different processes. Additionally, studies of resting state neuroelectric responses in the human brain have distinguished between modulations of oscillatory peaks from modulation of oscillatory troughs. That line of work has documented a difference between the variance of peaks and pits in restingstate time, (Amplitude Fluctuation Asymmetry Mazaheri10 ()), linking these to activity in visual areas, and suggesting this asymmetry derives from unbalanced forward vs. backward propagating currents within dendrites. Studies of functional Magnetic Resonance Imaging (fMRI) have also examined asymmetries in restingstate activity. These studies first identified all local minima and maxima in each time series and then contrasted the respective variances of the sets of local maxima and local minima. This variance asymmetry within spontaneous brain activity distinguishes adults from children Davis14 () and differentiates between wakefulness and sleep stages Davis16 (). As noted by Mazehari and Jensen Mazaheri10 (), “the amplitude fluctuations of oscillatory activity are conventionally viewed as being symmetric around zero”, but as we summarized above, emerging findings show that these conventions require revision, and importantly – that new and precise analytic tools are needed to quantify features of asymmetry within local minima and maxima. As we later show, amplitude asymmetry might be less than ideal in identifying such dynamics.
The need to capture and quantify possible asymmetries between the local maxima and minima of time series is not just inherent in neurobiological processes. Let us consider the time evolution of some financial index (e.g. S&P500) which represents the aggregate, collective evolution of the interaction of a number of financial assets over time. In quantitative finance, it is a wellknown empirical fact that qualitatively different dynamics operate microscopically when is on average increasing with respect to when this index is in a sustained decrease (which is equivalent to say that grows). In the former case, market is usually stable, pairwise correlation between the constituting assets is generally low (the system is said to be close to equilibrium Flanagan2016 ()) and risk perception is low, leading to a time series with low volatility. Conversely, a situation where decreases is indicative of a market under stress, where the correlation of the constituent assets grows due to common factors. As a result, any small and uncontrollable fluctuation can easily propagate throughout the system, hence the dynamics display larger volatility (larger uncertainty). The role that these two different market dynamics is playing can be therefore examined by analysing series statistics under index inversion , something that we will show is tightly related to asymmetries between peak and pit statistics.
Keeping in mind the aforementioned desiderata and as well as findings highlighting the importance of
discrete events, our aim here was to develop a general method for the efficient quantification of peak/pit asymmetry targeted at understanding the role
of local, nonoscillatory processes. After a thorough validation of such a methodology, our aim was then to apply it in a wide range of settings, including neurobiological, financial, climate time series and beyond. We capitalized on a recent general approach to the description of time series that provides information about both local and global temporal features, without assuming neither stationarity nor oscillations at any temporal scale.
This approach originates on the notion of a visibility graph Lacasa2008 (); Lacasa2009 (); multivariate (), a mapping that converts an ordered sequence of realvalued data (e.g. time series) into a graph of nodes, where each time point is mapped into a node and two nodes are connected if certain geometric and ordering criteria hold amongst the data (see figure 1 and Methods). Visibility Graphs were introduced with the aims of using the tools of Network Science Newmanbook () to describe the structure of time series and their underlying dynamics from a combinatorial perspective (other proposals for graphtheoretical time series analysis can be found in Kurths2017 () and references therein). Research on this methodology has been primarily mathematical (elaborating on mathematical methods severini (); Luque2016 (); nonlinearity (); motifs () to extract rigorous results on the properties of these graphs when associated to canonical models of complex dynamics epl (); jns (); pre2013 (); quasi (). In practice, this method can be used as a feature extraction procedure for constructing feature vectors for statistical learning purposes (see Shao2010 (); Ahmadlou2010 (); Ahmadlou2012 (); Bhaduri2016 (); Sanino (); meditation_VG () for just a few examples in the life sciences or physics3 (); fluiddyn0 (); fluiddyn1 (); fluiddyn2 (); physics2 (); suyal (); Zou (); physio1 () for other applications in the physical sciences).
Importantly, a fundamental property is that the conversion from time series to
visibility graphs is invariant under several transformations that map onto to common
nuisance (e.g., instrument) effects which are typical in neurobiological time series and beyond. These include linear trends,
amplitude modulation on longer scales, or variations (slowingdown, speedingup) in the rate of the process in question. Consequently, these invariances result in more efficient combination of information across
measurements after transforming a time series into a visibility graph. Visibility graphs provide informative features both at the local level of single
vertices and at the global level, where distributions of vertex features are described.
Such global features reflect, for instance, the selfsimilarity of fractal time series epl (), estimation of entropy production due to time irreversibility epjb (), discrimination between noise and chaos toral (); ravetti (), etc. Interestingly, previous works report that visibility graph features capture both linear and nonlinear information of the dynamical process and thus extend above and beyond standard power spectrumlike measures which only capture linear correlations.
Additionally, because visibility graphs can provide insights into local, nonoscillatory processes motifs (), they go
beyond the information provided by ‘global’ measures that summarize time
series in a single parameter NTSA () (series entropy, fractal dimensions, power spectrum
or even methods that are sensitive to similar (repetitive) motifs on multiple time scales (e.g., multiscale entropy methods)).
Equipped with the notion of visibility graph as a starting point, we explore here a systematic extension to that method, which is designed to satisfy the desiderata outlined above and capture different dynamics that may determine peaks and pits in a given signal. As indicated above, we are interested in a measure of asymmetry that can be applied to time series that do not have a natural minimum, which can be combined across measurements, and which is relatively robust to noise.
The rest of the paper proceeds as follows: in the Methods section we introduce the theoretical formalism, along with a theoretic analysis and validation for synthethic processes. For simple (stochastic and deterministic) processes we show that this methodology can capture peak/pit asymmetries with similar performance to that state of the art indicators, however for more complex processes involving combination of different dynamics and scales we show that this methodology outperforms state of the art approaches. We confirm such findings with extensive numerical simulations and rigorous results on concrete, canonical complex dynamics.
In the results section, we first show that this method offers novel descriptions for spontaneous brain activity in humans and differentiates between states of consciousness. For financial time series, it captures important features of financial regimes linked to major events in stock markets. We finally explore the application of this framework for extensive climatic data.
Ii Methods
ii.1 Graphtheoretical framework for peak/pit asymmetry quantification
Peak and Pit subsets.
Let us consider a realvalued time series of size , . The traditional step to get access to peak and pit statistics is to define two ordered sets, namely and .
The hypothesis is that the statistics of these sets, and its difference, carry relevant information on local fluctuations of and can be used as a feature for making diagnostics. Mathematically, differences in the statistics of peak and pit can be related in principle to two scenarios, namely: () different marginal distributions, and () different correlations (i.e. peak and pit might have similar marginals but different temporal correlations). Additionally, () characteristic crosscorrelations between peak and pit can also be informative (e.g. peak and pit might have similar marginal distributions but say, fluctuate in an anticorrelated way with one another). Among the plethora of possible descriptors, one should be able to identify and separate what scenarios () the measures are addressing. As a matter of fact, asymmetries between peak and pit statistics have only been addressed relatively recently, and the state of the art statistical tests e.g. Amplitude Fluctuation Asymmetry (AFA Mazaheri10 ()) or Amplitude Variance Asymmetry (AVA Davis14 ()) essentially considers scenario () by comparing the variances of these marginals (through the logarithm of the socalled Variance Ratio in the case of AVA, where is the variance of the random variable ).
As these measures only considers the onepoint marginals of each set, they cannot give us particular insights on scenarios () or (), and therefore disregard these aspects. However it is a computationally simple statistic, and probably because of its simplicity the quantity is currently used to assess the similarity between peak and pit statistics as mentioned above when describing prior applications to neurobiological time series. In addition, the degree to which AVA or AFA in any given time series is driven by extreme points Amor () is unclear and necessitates manual examination. In practice, it appears there may be a relation between the AVA of brain time series and the skewness of the empirical distribution Amor (), though analytically AVA and skewness are independent as the same distribution can produce time series with markedly different AVA values depending on specific sampling.
Incidentally, note that the quantity fulfils the axioms of a metric only in the case where is positive
Furthermore, note that by construction , meaning that this decomposition is lossy: in a generic fluctuating signal there are data which are neither in the peak nor in the pit set, so a priori important features of the local fluctuations might be lost if one only looks at the peak and pit sets and discard intermediate data.
An alternative, which we explore in what follows, is to be able to extract highorder features from peak and pit local neighborhoods by projecting the full signals into an appropriate topological space.
Visibility graphs and topbottom VG/HVG asymmetry (VGA). Consider again a time series . The socalled natural visibility graph (VG) Lacasa2008 () is extracted from the series by associating a different node to each datum, and linking every two nodes and with an edge if the following convexity criterion holds between the associated time series data:
Related to this graph, the socalled horizontal visibility graph (HVG) Lacasa2009 () is a subgraph of VG, obtained by applying a similar procedure with a slightly different linking criterion which only relies on the ordering of the data:
The theory of VG/HVG has been intensively used in recent years, not only to display a combinatorial representation of complex dynamics but also as a computationally efficient way of extracting informative topological features from empirical time series for statistical learning purposes. Among other properties, features of these graphs include a simple way to estimate the Hurst exponent in fractional Brownian motion epl () or the ability to easily distinguish fully chaotic processes from uncorrelated stochastic ones toral (); ravetti (); motifs () (both having identical flat power spectra). The Lyapunov exponent of a chaotic map can also be quantitatively obtained from the graph’s block entropies wolfram (). All in all, visibility graphs extracts valuable information from a time series on both linear and nonlinear level.
Here we label as (i.e. the ‘top’ visibility graph) to the graph extracted for either of the two aforementioned procedures
Our working hypothesis therefore exploits the potential differences between the top and bottom graphs as a proxy for quantifying the difference between local maxima and local minima statistics in , or the topbottom VG/HVG asymmetry (VGA) for short. Mathematically, VGA can be quantified in very many different ways (different protocols can be defined in an ad hoc way depending on the particular problem under study). For instance, simple global topological features that we know are informative include statistics of the degree sequence Luque2016 () such as graph’s degree distribution , for which, adopting the norm distance between distributions, leads to a particular definition of the asymmetry . This is just an informed choice and in any event, it always comes down to a comparison between certain set of features extracted from the top and bottom VG/HVG.
As mentioned above, there are a number of interesting scenarios involving differences in the marginal distributions or correlation structure of peaks and pits (, , ), and VGA can serve as a tool to investigate these asymmetries. While these issues have not been fully examined in prior work, one important study previous () has used a comparison between and with the intention of studying the different dynamics of local minima and maxima in the particular case of sunspot time series. As a matter of fact, sunspot series can be considered the degenerate case because they have a natural absolute and frequently occurring natural minimum (zero), which by definition imposes different features on the local minima (including impacting their variance, serial autocorrelation and power spectra). Similar cases where dampening function is applied () are similarly unhelpful.
Before presenting the practical algorithmic protocols we explore in synthethic processes the performance of this method, and in particular the ability to outperform current methods as well as the transversality (e.g. how well scenarios above are addressed).
ii.2 Theoretical validation on synthetic processes
To validate the method, we initially consider a battery of dynamical processes with varying degrees of complexity which focus on different aspects, and consider in every case the performance of a particular definition of (measured on both VG and HVG) and their comparison with the AVAbased metric (VR) (see figure 3 for an illustration). Results are summarized in table 1. We can highlight the following key results:

Symmetric stochastic processes (both with no correlations –white noise– and with rapidly decaying correlations –red noise– yield statistically identical top/bottom VG/HVGs and thus vanishing values of VGA.

For nonsymmetric white noise, peak and pit statistics are different, even if the process lacks information. In that case, both VGA (applied to VG) and VR detects such asymmetry, while VGA applied to HVG filters out dependences solely based on marginals and predicts no difference: in this sense we conclude that HVG does not capture scenario as defined above.
These two initial observations can be summarised in the following theorem:
Theorem 1. Let be a biinfinite time series generated by (i) white noise with a continuous probability density or by (ii) symmetric correlated red noise. Then VGA for HVG.
A proof for this theorem can be found in the appendix. 
For chaotic processes where peaks and pits are different both in terms of marginal statistics (scenario ) as well as in terms of correlation dependences (scenarios and ), all methods successfully detect such asymmetry. In particular, the following theorem holds:
Theorem 2. Let be a biinfinite time series generated by a fully chaotic logistic map. Then VGA for HVG.
A proof for this theorem can be found in the appendix. Gathering together theorems 1 and 2, we have rigorously proved that a chaotic process such as the fully chaotic logistic map can be easily distinguished from white noise, despite the fact that both processes have a flat power spectrum (deltadistributed autocorrelation function). Additionally, we numerically observe that chaotic processes with a fastdecaying correlation structure are also distinguishable from exponentially decaying correlated noise. 
Interestingly, there are several cases (such as processes with two alternating dynamics for peaks and pits) where currently used indicators (AVA, AFA) will fail, while VGA efficiently captures significant differences (i.e. VR does not capture scenarios and as defined above).

VGA is robust against noise pollution and works for short time series, enabling its use in practical cases.

VR is a quantity which by construction only depends on the marginal distribution of peaks and pits (in particular, the variance of these marginals), and as such it does not carry information on any temporal correlations, neither intra peaks or intra pits (scenario ), nor inter peaks/pits (scenario ). It is easy to prove that for a given series , if one breaks down all temporal correlations by reshuffling the series, then the new, reshuffled series and will still have the same marginal distributions and thus (under some assumptions) the same VR value, yet is indeed white noise with a flat spectrum and deltadistributed autocorrelation function, very different in general from the nonreshuffled series . In the same line, VR breaks down for any signal which is composed by two alternating processes with similar variances and dynamics with different invariance properties, whereas VGA is not afflicted by these drawbacks.
ii.3 Parametric analysis
Finally, in order to explore the relationship between VGA and a number of standard linear and nonlinear indicators (mean, variance, skewness, kurtosis, Lyapunov exponent), we have considered a parametric map (the logistic map) . As smoothly varies between and the map generates chaotic (aperiodic) time series with different degree of chaoticity (different Lyapunov exponent, attractor’s fractal dimension and physical invariant measure), interspersed by periodic windows. We recorded each statistic as we continuously scan .
In figure 4 we plot VGA (defined as the distance between top and bottom degree distributions for time series) extracted from the logistic map , for a range of values of for which the map generates chaotic trajectories of different type, interspersed by periodic windows. In the same figure we plot the Lyapunov exponent of the map. Both indicators are only weakly correlated (Pearson correlation coefficient ). Correlation between VGA and other standard indicators include: (mean), (standard deviation), (skewness), (kurtosis). We conclude that VGA finds negligible correlation with all indicators, except for standard deviation and Lyapunov exponent, where the correlation was found to be weak: VGA provides genuinely different information than standard linear and nonlinear indicators. Additionally, note that is notably unrelated to power spectrum (for instance, VGA is positive for and null for white noise, and both processes have identical power spectrum).
VGA  of :  
Process  Definition  Expected outcome  VG  HVG  AVA 
(symmetric) uniform white noise  ,  0  
exponentially correlated red noise  0  
Power law (nonsymmetric) white noise  ,  0  
Fully chaotic logistic map  
Fully chaotic logistic map polluted with white noise  , , , fully chaotic logistic map.  but decreasing with noise amplitude  
Selected pollution  , fully chaotic logistic map, ,  but decreasing with noise amplitude  
Alternating dynamics : fully chaotic logistic,  
Alternating dynamics : fully chaotic logistic,  where is such that the variance of the sinusoidal process is one  and independent of until both processes coalesce () 
From the theoretical analysis above, we conclude that VGA as defined above can be efficiently used to detect and quantify subtle differences between the statistics of local maxima and local minima in the associated series which are not correlated with standard linear and nonlinear indicators. If such asymmetry is solely based on the (onepoint) marginal distributions of peaks and pits (scenario ), then all three choices (VGA applied to VG and HVG, and AVA) are qualitatively similar. However, in the event the peaks and pits have identical marginals but different correlation structure, then only VGA capture this difference. Finally, if there is a difference between marginals but no difference in the correlation structure, then VGA (applied to VG) and AVA accurately capture the difference, but not VGA applied to HVG (as this latter is an order statistic).
Iii Data extraction methods and visibility protocol
We have applied the methodology in three different situations: EEGfMRI data, financial time series of NYSE, and worldwide temperature records. In what follows we provide details on data acquisition and preprocessing, as well as the protocols devised in each case to compute VGA.
iii.1 EEGfMRI data
Data acquisition and artifact correction
Data was collected during a synchronous EEGfMRI acquisition protocol, which was approved by the Ethical board of Goethe University (Kommission des Fachbereichs Medizin der J. W. GoetheUniversitat Frankfurt am Main as of January 10th, 2008). This data has been analyzed in several prior studies that examined: identification of Kcomplex correlates in N2 sleep Jahnke12 (), use of functional connectivity for sleep staging Tagliazucchi12 (), impact of sleep on serial autocorrelation of BOLD time series Tagliazucchi13 (), and our evaluation of AVA Davis16 (). Here we only used the EEG data to determine sleep stages. EEG was sampled at an initial sampling rate of 5 kHz, low pass filtered at 1 kHz, and downsampled to 250 Hz for artifact cleaning and sleep staging and further analysis (see below). Data were recorded using a BrainCapMR (Easycap) EEG cap (Herrsching, Germany) with 30 recording channels. The MRcompatible amplifiers were BrainAmp MR+, BrainAmp ExG; BrainProducts (Gilching, Germany). Data were recorded using BrainProduct’s “Recorder” and analyzed using Brain Product’s “Analyzer”. Analysis steps included: MRI and pulse artifact correction performed based on the Average Artifact Subtraction method Allen98 () as implemented in Vision Analyzer2 (Brain Products, Germany) followed by objective (CBC parameters, Vision Analyzer) ICAbased rejection of residual artifactladen components after artifact subtraction. Good quality EEG was obtained, which allowed for sleep staging by an expert, according to criteria of the American Academy of Sleep Medicine Berry12 (). Sleep staging was based on scoring 30s blocks of the EEG data. Based on this scoring we ignored sections associated with transitions between sleep stages, maintaining only those without transitions. BOLD time series sections matching these were spliced from the recorded time series. Functional MRI scans were acquired on a 3T system (Siemens Trio, Erlangen, Germany) using singleshot T2*weighted EPI (32 slices, repetition time/echo time , matrix , voxel size , distance factor ). To correct for physiological noise, physiological responses (cardiac, respiratory) were recorded through sensors from the MR scanner (sampling rate ) and MRcompatible devices (BrainAmp MR+, BrainAmp ExG; Brain Products). Sixtythree healthy nonsleepdeprived participants (thirtysix females, mean SD age of years) were scanned in the evening (starting from 8 PM). Data from those 55 participants who reached at least sleep stage N1 was used in our analysis.
fMRI epoch selection and preprocessing
Each of the 55 participants provided epochs of functional imaging data during wakefulness (W) and at least the N1 sleep stage. From these epochs, we set the minimal time series length at 135.2 seconds (65 volumes), amounting to 293 total epochs. In practice, most of the epochs were 510 times as long (see Davis16 () for details and quality control procedures).
We used AFNI Cox96 () for preprocessing and physiologicalnoise correction (PNcorrection). Despiking of the time series was carried out as part of the PNcorrection workflow using AFNI’s 3dDespike. This is important, because the procedure replaces time series that exceed 2.5Standard deviations of the mean with values estimated by smoothing to nearby values. This reduces the relative impact of extreme values.
The physiological (cardiac, respiratory) data were downsampled to the acquisition rate of single volume slices (15.4Hz). We used AFNI’s retroTS.m procedure to create slicebased regressors from these data. This utility generates, for each physiological regressor, a time series that is phase shifted to match the timing of each slice’s acquisition. From the cardiac and respiration recordings we derived 13 such slicebased regressors: 4 for the cardiac series and its harmonics, 4 for the respiratory series and its harmonics, and 5 for respiration variation over time and its harmonics, Birn06 (). We removed the variance explained by these 13 regressors from the BOLD time series using linear regression (RETROICOR; Glover00 ()) as implemented in AFNI.
Following PNcorrection, we discarded the first 5 volumes of each epoch from the analysis to allow for T1 stabilization effects, and then performed slice timing correction (3dTshift), motion correction (3dvolreg), and spatial smoothing (3dmerge, 6mm FWHM Gaussian Kernel) on all of the images. We then removed several sources of variance from the time series data via linear regression. These included (i) 6 motion parameters estimated during the head motion correction, and (ii) linear, secondorder and thirdorder polynomial trends. We used the residuals of the regression procedure as the functional MRI data of interest in all subsequent analyses.
iii.2 Calculating VGA and protocols
Neurobiological case
In a neurobiological context we were interested in monitoring systematic differences between regions of the degree distributions (rather than a net value out of the comparison of the whole distributions) which hold across participants in specific brain areas. A scalar projection such as the one defined for VGA in the synthetic cases lose such finegrained level of detail. Therefore, in this particular application we first constructed, for each voxel, degree distributions up to degree (this histogram truncation is justified as the modal degree value in fMRI series is around 3 or 4 with very few time points having a degree greater than 8, as shown in the results section). A second reason for selecting as a limit was our interest in local features of the time series, rather than in any highamplitude but infrequent spikes isolated from each other by more than 22 (TR = samples) seconds. We recall that Figure 1 illustrates graphically the concept of how node degrees are established, and how a time series can be characterized via degree distributions of the top and bottom VG/HVGs. As an additional preprocessing, in this application we created empirical cumulative distribution (CDF) functions of each histogram, normalizing the bar heights so that the area of the histogram is equal to one.
Normalization to common spatial reference space. After creating maps of visibility graph differences (topbottom) on a single voxel level, we obtained a transformation between each EPI time series to its corresponding anatomical image using FSL’s epireg script. The most important steps in this procedure are FAST’s (FMRIB’s Automated Segmentation Tool; Zhang01 () histogram based segmentation of the T1 structural scans to derive white matter maps, and the use of the boundaries of these white matter maps to perform BoundaryBased coRegistration of the EPIs to their corresponding T1 structural images (BBR; Greve09 ()). We then performed nonlinear normalization (FNIRT), of each subject’s T1 images into MNI space. Importantly, we concatenated the two transformations (EPI to T1; Linear, guided by white matter boundaries) and T1 to MNI (nonlinear warp) to derive a single transformation from EPI space to MNI space. We used this transformation to align the AVA maps from original space to the MNI template in a single step.
Singlevoxel calculations, and grouplevel analyses. The sleepstaging procedure allowed us to obtain fMRI time series for wakefulness (W) and three different sleep stages: N1, N2 and N3. On the single participant level, for each of the four conditions, we conducted voxelwise analyses to derive the nodedegree histograms for the top and bottom graphs (i.e. these were derived for each voxel). These were represented as the empirical cumulative distribution function. With these CDFs we could answer the following two questions:

On the singlecondition level we identified brain areas that differentiated the top from bottom CDFs; we refer to this ‘difference’ histogram as Asymmetry of Visibility Histogram (AoVH), which is our primary way of implementing VGA in this application. We performed these analyses for CDFs derived from both VG and HVG. It is important to note that because these CDFs were normalized, a main effect of condition is not possible (the histograms average across bins to one) and there will always be a strong main effect of bins. Instead, it is the interaction between top/bottom analysis and Bin which we use to show differences across conditions. In reference to the illustration in Figure LABEL:fig:VHGVG (panel C), we are interested in whether the difference between the blue and red histograms, calculated for a given voxel’s time series’ are systematic across participants.

Then, in order to study differences between conditions, we compared the AoVH of the two conditions. We did this by deriving AoVH for each condition and then determining statistically whether these differed between the two conditions. Note that comparisons between conditions were based on time series matched for length within participant.
On the group level, to evaluate voxels showing statisticallysignificant AoVH within each study condition (W, N1, N2, N3) we conducted voxelwise repeated measures ANOVA with two fixed factors; Time series view (two levels [Top, Bot]), and nodedegreehistogram Bin (8 levels as explained above). To enable inferences about the population, we modeled participants as a random factor. In order to compare between any two conditions (e.g., W vs. N1) we derived the AoVH for each condition, and conducted a repeated measures ANOVA with 2 fixed factors: Condition (e.g., W vs. N1) and Node degree histogram bin of the difference histogram (8 bins, values 2(min)9(max)). We note that 8 bins are used here, as the 9th bin would not be independent of the prior 8 (due to normalization).
These analyses returned a statistical significance value for the interaction term, for each brain voxel. We then implemented a clustering procedure forman99 () to identify brain areas where many contiguous voxels, each with a value of .001 show a significant interaction: this identifies an ‘activation cluster’.
Financial case
In application to financial data, we considered a dataset of financial stocks comprising stock evolution of 35 major American companies from the New York Stock Exchange (NYSE) and NASDAQ in the period 19982012, the majority of which belong to the Dow Jones Industrial Average. Data have been extracted from Flanagan2016 ().
In previous sections we have confirmed that the () distance between top and bottom degree distributions of VG/HVG is an efficient and informative definition of VGA which in many cases outperforms state of the art signal peak/pit asymmetries and is genuinely different from standard linear and nonlinear indicators. This is precisely the scalar which we shall use in the financial context. Our protocol is as follows: we compute the VGA by splitting the time series for each company into yearly time series and calculate the distance over the associated HVGs for each year. Accordingly, each year is characterized by a vector of distances. For instance, for year 1998, we have a vector of 35 dimensions whose entry is the VGA for company in 1998. Subsequently, a principal component analysis is performed to dimensionally reduce data, and a projection into the two first principal components is used to visually cluster different years.
Climatic case
This application is conceptually similar to that reported for financial time series. Rather than considering stock prices of major shares in multiple years, we considered daily temperature data sampled over a global grid with a resolution of 192 longitude and 94 latitude, with average temperature for each day (365 values) provided for each grid point (data obtained from Kanamitsu2002 ()). For each year between 1995 and 2015 we compute VGA following the same definition used in the financial setting.
Iv Results
iv.1 Application to fMRI
We first compared the top and bottom visibility graphs within the Wakefulness and the three sleep stages (N1, N2 and N3). As shown in Figure 5, we found significant differences, predominantly in thalamic and frontal regions, for each of these conditions, with the exception of the case of VG in the N3 sleep condition.
To better understand these results, we determined which visibility values tended to contribute most strongly to the statisticallysignificant differences in degree distributions that produced the Wakefulness (W) results in Figure 5 (rows 1, 2). To this end, for each cluster we derived a histogram that communicated the visibility bins that most strongly differentiated the top and down degree distributions for each voxel in the cluster. We did this by (i) transforming the cluster to original space, (ii) for each voxel, identifying the bin that maximally differentiated the top from down histogram, (iii) transforming that value to common space, (iv) creating an average across participants for each voxel in the cluster. The resulting histograms communicated a very clear and consistent result: for VG the modal degree value that maximally differentiated the top and bottom histograms was 4, with narrow tails towards the values of 3 and 5 (indicating that for some participants, some voxels maximally differentiated the histograms at values 3 and 5). Importantly, there were no cases with means below 3 or above 5. For HVGs, in all clusters the modal degree value that maximally differentiated the top and bottom histograms was 3, with very narrow tails towards 2 and 4.
These findings are very important as they show that the differences identified by VGA (AoVH) were driven by very local dynamics rather than due to differences in propensity of isolated extreme events.
In a separate analysis, we also found that VGA profiles could discriminate wakefulness from sleep. As shown in Fig.6, we identified numerous areas, in both occipital (visual cortex) and lateral temporal cortices, where dynamics during wakefulness and N2 sleep differed significantly. No differences were found between W and N1 or between W and N3.
iv.2 Application: unsupervised clustering of financial periods
First, we show in the upper panels of Fig.7 that VGA is not correlated with the measure based on Variance Ratio (VR, left panel, Pearson correlation ), nor with the standard econometric measure to capture financial series fluctuations (annualized volatility, right panel, Pearson correlation ). This means that VGA is a genuinely different measure in this domain. The correlation between an index and its volatility is indeed related to the skewness of the index, hence this result further validates the hypothesis that VGA is independent of skewness.
Then, we perform principal component analysis (PCA) on the vector where is the VGA (HVG) for company for a specified year, and project results on a twodimensional space spanned by the first two principal vectors of the PCA projection. This plot is shown in the left panel of Figure 7,
and indicates that the measure based on VGA is informative and we can cluster periods of financial turmoil together (the global financial crisis 20092010 is found separated from the dotcom bubble 19982000 and from the rest of the years).
In the right panel of the same figure, we further compare equivalent plots produced via . We thus compute for each year, for each company, and perform PCA on the space . In this case no robust clustering appears.
The interpretation of these results is as follows: the local fluctuations in the time arrangement of local maxima differs from the arrangement of local minima, and these differences vary with respect to the overall state of financial stress, i.e. collectively the set of stock prices is responsive to the financial stress state, and the difference between peak and pit statistics is informative of such collective financial state. As a result, it is easy to cluster apart years where the financial system was ’in equilibrium’ from periods where the financial system is under stress and driven out of equilibrium, such as periods of financial bubbles or global crisis, which emerge as different periods in PCA space. We find that the statistic fails to capture such traits: this is, according to our previous theoretical analysis, indicative that the difference between local maxima and local minima which cluster apart financial periods are not simply related to differences between marginal distributions, but rely on more subtle differences in the correlation (temporal ordering) structure.
iv.3 Application: global daily temperature time series
The VGA map for the year 2015 is seen in Panel A of 8. As a first step we quantified the relation between VGA and several moments of the temperature distribution. To this end, we derived VGA spatial heat maps for each of the years 19952015, and we then calculated the correlation between the observed value of VGA (per grid point) and the following parameters of the yearly temperature time series: (i) mean, (ii) standard deviation, (iii) skewness , and (iv) kurtosis. Note that this returns a set of pairwise correlation values per year and we can then examine the distribution of these correlation values across years. Similarly to what we found for the validation case, the correlation with mean was low (Mean Pearson’s R; ; the correlation with standard deviation, r = 0.24 ; with skewness, r = 0.09, and with kurtosis, r = 0.13, . Note the low correlation with skewness, which strongly suggests that VGA measure is not loading on extreme events Amor (). The correlation between VGA and the standard deviation was moderate, and the correlations with mean and kurtosis minimal. This suggests that VGA provides information not captured by typical moments.
A Principal Component Analysis applied to the 21 Distance maps (19952015) identified a first component that accounted for 40% of the variance, and a second that accounted for 5% (see 8). The yearly loadings on the first component were quite stable across the years. In contrast, the yearly loadings on the second component show a strong linear change with time.
V Discussion
In this work, we first presented the theoretical and practical motivations for developing new methods for studying detailed features of local minima and maxima in various temporal domains. We provide analytical arguments that suggest that such issue can be addressed by a combinatorial method –originally studied in the context of solar activity previous ()– based on deriving two visibility graphs from a single time series. We have thoroughly validated the concept by proposing a detailed methodology which we apply on a battery of synthetic time series extracted from complex dynamics of different origin (encompassing correlated stochastic processes, chaotic dynamics and processes involving multiple scales). We have shown, both numerically and rigorously, that exploiting the difference and asymmetry between top and bottom VG/HVG (what we labeled as VGA) correctly identifies the peak/pit asymmetries and outperforms standard methods in several cases.
In applying this method, we first found that it offered a new view into spontaneous resting state dynamics in the human brain. While prior neuroimaging work based on amplitudevarianceasymmetry had pointed to sensory cortices as ones having different peak vs. pit dynamics AVA (), the current results identify frontal regions as exhibiting asymmetric restingstate BOLD fluctuations during wakeful rest. Furthermore, these asymmetric patterns were also found during the deeper sleep stages (N2 and N3), which is a departure from prior findings (DAVIS16) where AVA failed to identify such signatures. An important result was that in all clusters, these dynamics were driven by differences in the relative frequency of time points (nodes) with relatively low degree, indicating that these differences are not due to a difference in relative frequency of rare, extreme events but due to differences in frequency of time points with relatively moderate connectivity – i.e., a very local phenomenon.
Then, we were able to cluster periods of financial activity according to differences in peak/pit asymmetry via principal component analysis in visibility graph feature space. Clusters gathering periods of similar financial stress emerge naturally, and capture finergrained structure than a basic analysis based on AVA. Finally, we demonstrated a straightforward application of the methods in question to daily temperature time series. The latter produced two main results. First, already hinted at by our analysis of synthetic series, there were either very modest or no correlations between VGA and typically studied distribution moments including mean, standard deviation, skewness or kurtosis. Second, A PCA analysis revealed well structured spatiotemporal correlations between changes in VGA in several latitude bands, broadly matching planetaryscale topology of the subtropical and polar jet streams. The yearly PCA loadings of the main component did not identify a strong cyclical pattern consistent with ElNino events, but the loadings on component 2 suggest a gradual monotonic change in VGA correlations over time between several areas.
In addition, we have shown that VGA goes beyond basic statistics defined over VG/HVG since topological features of a standard VG/HVG cannot capture accurately the situation when the signal is aperiodic and alternates between large and small values (i.e. the structure of small data is completely lost as they will always be nodes with degree ), such being for instance the case of a chaotic process with a disconnected twoband chaotic attractor among other cases.
We anticipate that VGA will provide a efficient feature extraction method which will be valuable for statistical learning analysis across the disciplines. Finally, note that in this work the definitions of VGA mainly reduce to a comparison between global topological features. However the same methodology also allows to extract local topological features (e.g. motifs motifs ()), whose comparison can thus be exploited for early warning and anomaly detection purposes.
Vi APPENDIX: Proofs of the theorems
Theorem 1. Let be a biinfinite time series generated by (i) white noise with a continuous probability density or by (ii) symmetric correlated red noise. Then VGA for HVG.
Proof. The proof for this theorem is straightforward. Firstly, we consider the case (i). We rely on a result of Lacasa2009 () where it has been proved that the degree distribution of the HVG of generated by white noise is
(1) 
Trivially, if is white noise extracted from , this means that it is a collection of i.i.d. random variables extracted from . But then is again a collection of i.i.d. variables extracted from a different density . By virtue of a theorem proved in ?, the degree distribution of the HVG associated to such series coincides with eq.1, hence VGA. This closes the first part of the proof.
For (ii), we only need to remark that red noise (AR(1)) is generated via the stochastic map
Trivially, if is drawn from a symmetric distribution (such as the Gaussian function above), then both
and are equally likely realizations of eq. VI, hence VGA vanishes. ∎
Theorem 2. Let be a biinfinite time series generated by a fully chaotic logistic map. Then VGA for HVG.
Proof. The proof for this theorem is somewhat more convoluted than for theorem 1. Since VGA is semipositive definite, it vanishes if and only if . Thus, in order to prove positivity we will only need to find that for a given degree , . After a quick numerical inspection, we choose . In toral () it was analytically proved that for a fully chaotic logistic map, . Without loss of generality, we will prove that .
We capitalise on the perturbative expansion formalism advanced in nonlinearity () to analytically compute the degree probabilities of HVG associated to chaotic maps with smooth physical invariant measure. This expansion is diagrammatic, and in particular it expands on the number of hidden nodes (see nonlinearity () for details). In order to apply this technique here, we first observe that the bottom HVG of is equivalent to the top HVG of . Accordingly, the set of diagrams that comply with in the original series is schematized in figure 9. It is easy to see that this is an infinitely unfolding combination of diagrams (each with a different number of hidden nodes), and as such formally we have
where is the number of hidden nodes in each diagram. Let us denote as the fully chaotic logistic map, as the invariant measure of the map. Since this is a Markovian and deterministic map, the probability of observing a datum given that the previously observed datum is is simply , where is the Diracdelta distribution. The integral corresponding to each contribution to the probability is then:
(2)  
The first important point to note is that the first integral can be computed using the scaling properties of the Dirac delta:
where are the roots of the equation (in particular, we only need to sum up over those roots that belong to the interval of integration, that is, only those roots that satisfy . For , the only root which fulfils such inequality is , hence
(3) 
The rest of the integrals only have the effect of reducing the interval of integration of , as these are Dirac integrals and thus are either one or zero, since is one if and zero otherwise.
Let us assume for a moment that we remove in the integral above any contribution of nodes above . After a little algebra, this particular case labeled would yield a probability
(4) 
where was defined in eq.3, and the shrinking of the interval of integration in from to comes from the integral , which is only equal to one if , what happens for (see figure 10). Importantly, any other additional Dirac delta integral (associated to the contribution of further nodes –hidden or otherwise–) of the form will have the effect of further shrinking the interval in the integral over , in such a way that for each order the original interval of integration of shrinks into a subinterval . To showcase such situation let us consider in full detail the computation of and .
Let’s consider the last two integral terms above: they require
and . Looking at figure 10 and considering the fixed point structure of and higher iterates of the map, we find that the first condition shrinks the interval of integration of from to , whereas the second condition is fulfilled for . The intersection of the two conditions above yields an interval
, and therefore
Let’s then consider the last three integral terms above: they require
, , . Looking at figure 10 and considering again the fixed point structure of and higher iterates of the map, we find that the first condition shrinks the interval of integration of from to , whereas the second condition is fulfilled now for
. The third condition is fulfilled in
. The intersection of all these conditions thus provides
, and therefore
One can proceed indefinitely adding higher orders in , which yield smaller and smaller contributions to the total probability. In other words, for any the resulting interval is always a subset of . Since the integrand in the last integral of eq.4 is always positive (see figure 11 for a graphical support), we conclude that . Furthermore, we need to prove that if we concatenate the subinterval obtained for each , its union will always be smaller than :
(5) 
This last statement can be easily proved by showing that there are always subintervals of not covered in any . In particular, consider the subinterval . We have explicitely shown that such subinterval is not included in , nor in . Now, for by construction one of the conditions which is always present is . Now, such condition is fulfilled for , leaving our interval out. As this situation holds , this directly yields eq.5.Finally, since integration is monotonic and since the integrand is positive, we trivially have
Altogether, this means that , what concludes the proof. ∎
Acknowledgments. We wish to thank Dr. Alon Angert for his advice on meteorological time series. UH‘s work was conducted in part while serving at and with support of the National Science Foundation. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the NSF. LL’s acknowledges funding from an EPSRC Early Career Fellowship EP/P01660X/1.
Author contributions. UH and LL designed research; UH, BD, JI, RF and LL conducted research, ET and HL contributed data, all authors discussed results and UH, JI, RF and LL wrote the paper.
Competing financial interests. None.
Footnotes
 This can be seen from the triangle inequality: take , , . For we can drop the absolute values and the triangle inequality saturates Take however , . In that case the triange inequality is not satisfied in general, as we have which in general is not larger or equal than .
 at this point we consider both kinds of graphs VG and HVG indistinctively, although we know that they indeed provide different sorts of information and in practice depending on the particular processes under study it will be more adequate to make use of either HVG or VG
 For another simple example, consider two time series, the first being a periodic series of period 2, and the second being a mix of two processes: uniform white noise in for even times , and a constant value for odd times . The associated VG/HVG is identical despite that both processes are qualitatively different, however using both top and bottom graphs the difference is obvious.
References
 Kantz, H., & Schreiber, T. Nonlinear time series analysis (Vol. 7, Cambridge university press, 2004).
 Tsay, R. S. Analysis of financial time series (Vol. 543, John Wiley & Sons, 2005).
 L. Lacasa, B. Luque, F. Ballesteros, J. Luque, J.C. Nuño, From time series to complex networks: the visibility graph, Proc. Nat. Acad. Sci. USA 105(13) (2008) 49724975.
 B. Luque, L. Lacasa, F. Ballesteros, J. Luque, Horizontal visibility graphs: Exact results for random time series, Phys. Rev. E 80(4) (2009) 046103.
 L. Lacasa, V. Nicosia, V. Latora, Network Structure of Multivariate Time Series, Sci. Rep. 5, 15508 (2015)
 M. Newman, The structure and function of complex networks, SIAM Review 45, 167256 (2003).
 L. Lacasa, B. Luque, J. Luque and J.C. Nuño, The Visibility Graph: a new method for estimating the Hurst exponent of fractional Brownian motion, EPL 86 (2009) 30001.
 B. Luque, L. Lacasa, Canonical horizontal visibility graphs are uniquely determined by their degree sequence, Eur. Phys. J. Sp. Top. 226, 383 (2017).
 L. Lacasa, On the degree distribution of horizontal visibility graphs associated to Markov processes and dynamical systems: diagrammatic and variational approaches, Nonlinearity 27 (2014) 20632093.
 S. Severini, G. Gutin, T. Mansour, A characterization of horizontal visibility graphs and combinatorics on words, Physica A 390, 12 (2011) 24212428.
 L. Lacasa and w. Just, Visibility graphs and symbolic dynamics, arXiv:1704.06467
 Review by Kurths on EPL 2017
 L. Lacasa and R. Toral, Description of stochastic and chaotic series using visibility graphs, Phys. Rev. E 82, 036120 (2010).
 Ravetti, M. G., Carpi, L. C., GonÃ§alves, B. A., Frery, A. C., & Rosso, O. A., Distinguishing noise from chaos: objective versus subjective criteria using horizontal visibility graph. PLoS One, 9, 9 (2014), e108004.
 B. Luque, L. Lacasa, F. Ballesteros, A. Robledo, Analytical properties of horizontal visibility graphs in the Feigenbaum scenario, Chaos 22, 1 (2012) 013109.
 B. Luque, A. Núñez, F. Ballesteros, A. Robledo, Quasiperiodic Graphs: Structural Design, Scaling and Entropic Properties, Journal of Nonlinear Science 23, 2, (2012) 335342.
 A.M. Núñez, B. Luque, L. Lacasa, J.P. Gómez, A. Robledo, Horizontal Visibility graphs generated by typeI intermittency, Phys. Rev. E, 87 (2013) 052801.
 Lacasa, L., Nunez, A., Roldan, E., Parrondo, J. M., & Luque, B, Time series irreversibility: a visibility graph approach. Eur. Phys. J. B, 85, 6 (2012) 217.
 Y Zou, M. Small, Z. Liu and J. Kurths, Complex network approach to characterize the statistical features of the sunspot series, New J. Phys. 16 (2014).
 Shao, Z. G. Network analysis of human heartbeat dynamics. Applied Physics Letters, 96, 7 (2010) 073703.
 Ahmadlou M, Adeli H, Adeli A. New diagnostic EEG markers of the Alzheimer’s disease using visibility graph. Journal of neural transmission 117(9) (2010)1099109.
 Ahmadlou, M., Adeli, H., & Adeli, A. Improved visibility graph fractality with application for the diagnosis of autism spectrum disorder. Physica A 391, 20 (2012) 47204726.
 Bhaduri, S., & Ghosh, D. Electroencephalographic data analysis with visibility graph technique for quantitative assessment of brain dysfunction. Clinical EEG and neuroscience 46, 3 (2015), 218223.
 Sannino, S., Stramaglia, S., Lacasa, L., & Marinazzo, D. Visibility graphs for fMRI data: multiplex temporal graphs and their modulations across resting state networks. Network Neuroscience (2017)
 S. Jiang, C. Bian, X. Ning and Q.D.Y. Ma, Visibility graph analysis on heartbeat dynamics of meditation training, Appl. Phys. Lett. 102 253702 (2013).
 Davis, B., Jovicich, J., Iacovella, V., & Hasson, U. (2013). Functional and developmental significance of amplitude variance asymmetry in the BOLD restingstate signal. Cerebral Cortex, 24(5), 13321350.
 Allen, P. J., Polizzi, G., Krakow, K., Fish, D. R., & Lemieux, L. Identification of EEG events in the MR scanner: the problem of pulse artifact and a method for its subtraction. NeuroImage, 8(3) (1998), 229239.
 Berry, R. B., Brooks, R., Gamaldo, C. E., Harding, S. M., Marcus, C., & Vaughn, B. The AASM manual for the scoring of sleep and associated events. Rules, Terminology and Technical Specifications, Darien, Illinois, American Academy of Sleep Medicine. (2012).
 Cordes, D., Haughton, V. M., Arfanakis, K., Carew, J. D., Turski, P. A., Moritz, C. H., .. & Meyerand, M. E. Frequencies contributing to functional connectivity in the cerebral cortex in ”restingstate” data. AJNR Am J Neuroradiol, 22(7), 2001, 13261333.
 Cox, R. W. AFNI: software for analysis and visualization of functional magnetic resonance neuroimages. Comput Biomed Res, 29(3), 1996, 162173.
 Davis, B., Jovicich, J., Iacovella, V., & Hasson, U. Functional and developmental significance of amplitude variance asymmetry in the BOLD restingstate signal. Cerebral Cortex, 24(5), 2014, 13321350. doi: 10.1093/cercor/bhs416
 Davis, B., Tagliazucchi, E., Jovicich, J., Laufs, H., & Hasson, U. Progression to deep sleep is characterized by changes to BOLD dynamics in sensory cortices. NeuroImage, 130, 2016, 293305. doi: 10.1016/j.neuroimage.2015.12.034
 Fischl, B., Sereno, M. I., Tootell, R. B., & Dale, A. M. Highresolution intersubject averaging and a coordinate system for the cortical surface. Hum Brain Mapp, 8(4), 1999, 272284.
 Glover, G. H., Li, T. Q., & Ress, D. Imagebased method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn Reson Med, 44(1), 1999, 162167.
 Greve, D. N., & Fischl, B. Accurate and robust brain image alignment using boundarybased registration. NeuroImage, 48(1), 2009, 6372. doi: 10.1016/j.neuroimage.2009.06.060
 Jahnke, K., von Wegner, F., Morzelewski, A., Borisov, S., Maischein, M., Steinmetz, H., & Laufs, H. To wake or not to wake? The twosided nature of the human Kcomplex. NeuroImage, 59(2), 2012, 16311638. doi: 10.1016/j.neuroimage.2011.09.013
 Liu, X., & Duyn, J. H. Timevarying functional network information extracted from brief instances of spontaneous brain activity. Proc Natl Acad Sci U S A, 110(11), 2012, 43924397. doi: 10.1073/pnas.1216856110
 Mazaheri, A., & Jensen, O. Rhythmic pulsing: linking ongoing brain activity with evoked responses. Front Hum Neurosci, 4, 2010, 177, doi: 10.3389/fnhum.2010.00177
 Patel, A. D., & Balaban, E. Temporal patterns of human cortical activity reflect tone sequence structure. Nature, 404(6773), 2000, 8084. doi: 10.1038/35003577
 Petridou, N., Gaudes, C. C., Dryden, I. L., Francis, S. T., & Gowland, P. A. Periods of rest in fMRI contain individual spontaneous events which are related to slowly fluctuating spontaneous activity. Hum Brain Mapp, 2013, 34(6), 13191329. doi: 10.1002/hbm.21513
 Tagliazucchi, E., Siniatchkin, M., Laufs, H., & Chialvo, D. R. The VoxelWise Functional Connectome Can Be Efficiently Derived from Coactivations in a Sparse SpatioTemporal PointProcess. Front Neurosci, 2016, 10, 381. doi: 10.3389/fnins.2016.00381
 Tagliazucchi, E., von Wegner, F., Morzelewski, A., Borisov, S., Jahnke, K., & Laufs, H. Automatic sleep staging using fMRI functional connectivity data. NeuroImage, 63(1), 2012, 6372. doi: 10.1016/j.neuroimage.2012.06.036
 Tagliazucchi, E., von Wegner, F., Morzelewski, A., Brodbeck, V., Jahnke, K., & Laufs, H. Breakdown of longrange temporal dependence in default mode and attention networks during deep sleep. Proceedings of the National Academy of Sciences, 2013, 110(38), 1541915424.
 Zhang, Y., Brady, M., & Smith, S. Segmentation of brain MR images through a hidden Markov random field model and the expectationmaximization algorithm. IEEE transactions on medical imaging, 2001, 20(1), 4557.
 Birn, R. M., Diamond, J. B., Smith, M. A., & Bandettini, P. A. Separating respiratoryvariationrelated fluctuations from neuronalactivityrelated fluctuations in fMRI. NeuroImage, 2006, 31(4), 15361548.
 Forman, S.D., Cohen, J.D., Fitzgerald, M., Eddy, W.F., Mintun, M.A., Noll, D.C. Improved assessment of significant activation in functional magnetic resonance imaging (fMRI): use of a clustersize threshold. itMagn Reson Med, 1995, 33(5), 63647.
 A. Aragoneses, L. Carpi, N. Tarasov, D.V. Churkin, M.C. Torrent, C. Masoller, and S.K. Turitsyn, Unveiling Temporal Correlations Characteristic of a Phase Transition in the Output Intensity of a Fiber Laser, Phys. Rev. Lett. 116, 033902 (2016).
 M. Murugesana and R.I. Sujitha1, Combustion noise is scalefree: transition from scalefree to order at the onset of thermoacoustic instability, J. Fluid Mech. 772 (2015).
 A. Charakopoulos, T.E. Karakasidis, P.N. Papanicolaou and A. Liakopoulos, The application of complex network time series analysis in turbulent heated jets, Chaos 24, 024408 (2014).
 P. Manshour, M.R. Rahimi Tabar and J. Peinche, Fully developed turbulence in the view of horizontal visibility graphs, J. Stat. Mech. (2015) P08031.
 RV Donner, JF Donges, Visibility graph analysis of geophysical time series: Potentials and possible pitfalls, Acta Geophysica 60, 3 (2012).
 V. Suyal, A. Prasad, H.P. Singh, VisibilityGraph Analysis of the Solar Wind Velocity, Solar Physics 289, 379389 (2014)
 Y. Zou, R.V. Donner, N. Marwan, M. Small, and J. Kurths, Longterm changes in the northsouth asymmetry of solar activity: a nonlinear dynamics characterization using visibility graphs, Nonlin. Processes Geophys. 21, 11131126 (2014).
 J.F. Donges, R.V. Donner and J. Kurths, Testing time series irreversibility using complex network methods, EPL 102, 10004 (2013).
 R. Flanagan and L. Lacasa, Irreversibility of financial time series: a graphtheoretical approach, Physics Letters A 380, 16891697 (2016)
 M. Kanamitsu, W. Ebisuzaki, J. Woollen, SK Yang, J.J. Hnilo, M. Fiorino, G. L. Potter. NCEPDOE AMIPII Reanalysis (R2) Bulletin of the American Meteorological Society 16311643. 2002.
 J. Iacovacci and L. Lacasa, Sequential visibilitygraph motifs, Phys. Rev. E 93, 042309 (2016)
 Amor, T. A., Russo, R., Diez, I., Bharath, P., Zirovich, M., Stramaglia, S., Cortes J.M., De Arcangelis L., & Chialvo, D. R. Extreme brain events: Higherorder statistics of brain resting activity and its relation with structural connectivity. EPL (Europhysics Letters) 111, 6 (2015), 68007.