Assessment of visibility graph similarity as a synchronization measure for chaotic, noisy and stochastic time series

Finding synchronization between the outputs of a dynamic system, which are represented mostly as time series, helps to characterize the system activities during an occurrence. An important issue in analyzing time series is that they may behave chaotically or stochastically. Therefore, applying a reliable synchronization measure which can capture the dynamic features of the system helps to quantify the interdependencies between time series, correctly. In this paper, we employ similarity measures based on visibility graph (VG) algorithms as an alternative and radically different way to measure the synchronization between time series. We assess the performance of VG-based similarity measures on chaotic, noisy and stochastic time series. In our experiments, we use the Rössler system and the noisy Hénon map as representative instances of chaotic systems, and the Kuramoto model for studying detection of synchronization between stochastic time series. Our study suggests that the similarity measure based on the horizontal VG algorithm should be favored to other measures for detecting synchronization between chaotic and stochastic time series.


Introduction
In recent decades, various local and global techniques, operating in time, frequency, or wavelet domain, have been introduced for measuring synchronization among time series. Local techniques establish synchrony between pairs of time series, whereas global techniques can handle an arbitrary number of time series (Dauwels et al. 2010). The crosscorrelation function and its counterpart in the frequency domain, i.e., the coherence function, were the first linear methods developed for quantifying synchronization between time series (Barlow and Brazier 1954;Brazier 1968). These were followed by the development of nonlinear techniques based on mutual information (Panzeri et al. 1999), nonlinear regression (da Silva et al. 1989), and phase synchronization (Rosenblum et al. 2004) among others, summarized overviews of which can be found in several reviews (Dauwels et al. 2010;Pereda et al. 2005;Ansari-Asl et al. 2006).
Recently, the visibility graph (VG) algorithm has attracted attention as a capable technique for time series analysis (Lacasa et al. 2008;Lacasa and Toral 2010). A time series (or signal) can be considered as an ordered set of values. By applying a mapping technique, we can turn this set into a different mathematical object (e.g., a graph or network) and investigate which properties of the time series are conserved, transformed or missed. Therefore, the preceding mapping gains some practical interest as it opens the possibility of analyzing a time series from an alternative angle. The motivation is completed when the mapping technique belongs to a relatively mature mathematical field, where information encoded in such a representation can be effectively extracted and processed.
Furthermore, time series analysis has some limitations, when it refers to study nonlinear or complex time series. Dynamical phenomena such as chaos, long-range correlated stochastic processes, intermittency and multifractality are examples of complex phenomena where time series analysis is pushed to its own limits. A new approach, like VG algorithm, that deal with the intrinsic nonlinearity by being intrinsically nonlinear are not only welcome, but also needed. Moreover, the technological era brings us the possibility of digitally analyzing massive data sets. With the aid of pf a well-suited algorithm we can filter such plethora of data. The family of visibility algorithms constitute one of other possibilities to map a time series into a graph and subsequently analyze the structure of the series through the set of tools developed in the graph / complex network theory (Lacasa et al. 2008;Lacasa and Toral 2010;Luque et al. 2009;Nuñez et al. 2012).
The idea of mapping time series into graphs seems attractive because it lays a bridge between nonlinear signal analysis and complex networks theory. Hence, the visibility graph (VG) algorithm has attracted attention as a capable technique for time series analysis due to its intrinsic non-locality, low computational cost, simple rules and straightforward implementation (Lacasa et al. 2008;Lacasa and Toral 2010;Luque et al. 2009;Nuñez et al. 2012). The VG algorithm provides an effective method to map time series to a graph permitting a mutual relationship between dynamical properties of time series and topological properties of the graph. Therefore, the information on time series is obtained just by analyzing the characteristics of the graph. In particular, it has been shown that both the structure of complex, irregular time series and nontrivial ingredients of its underlying dynamics are inherited in the topology of the visibility graphs, and therefore simple topological properties of the graphs can be used as time series features for description and automatic classification purposes. Examples include a topological characterization of chaotic series and the method has been used for the description and classification of empirical time series appearing in physics, physiology, neuroscience or finance (Lacasa et al. 2008;Lacasa and Toral 2010;Luque et al. 2009;Nuñez et al. 2012;Czechowski et al. 2016;Shao 2010;Zhu et al. 2012;Ahmadi and Pechenizkiy 2016;Yang et al. 2009).
For synchronization purpose, the series of connectivity degree (i.e., the number of edges connected to a node) of the visibility graph nodes is considered as a new time series, which is called the degree sequence (DS) time series. Measuring the synchronization between two DSs is called the VG similarity and can be presented as an alternative measure of synchronization between time series (Ahmadi and Pechenizkiy 2016;. An important issue in analyzing a dynamic system is that the system may behave chaotically and/or stochastically (Lacasa and Toral 2010;Korn and Faure 2003). Chaotic systems display sensitivity to an initial condition which manifests instability everywhere in the phase space and leads to non-periodic time series. When one or more parts of a dynamic system have randomness, it is called a stochastic system. A stochastic system does not always produce the same output for a given input.
In this paper, we present the main results of our experimental study aimed at assessing the capability of the VGbased similarity techniques in measuring synchronization between chaotic, noisy and stochastic time series. For this purpose, we conducted a comprehensive evaluation of the VG-based similarity measures on coupled Rössler system, noisy Hénon map, and the Kuramoto model. We compare the performance of VG-based similarity techniques with other commonly used synchronization measures, including cross-correlation, synchronization likelihood and variants of coherence and phase lag index. The results of our study suggest to choose the horizontal VG-based similarity measure for detecting synchronization between chaotic and stochastic time series.
The rest of this paper is organized as follows: we introduce the VG-based similarity measures in Sec. 2, discuss the experimental results in Sec. 3 and conclude with final remarks in Sec. 4.

Visibility graph algorithm
Let y i be a univariate time series of N data (i = 1,2,…,N). The VG algorithm converts the time series y i to a graph G(v i ), as a data point y i is mapped into a node v i in G. The signal node (i.e., t i , a point on the time series) represents a moment in which the data is recorded. Therefore, a time series of size N maps to a graph with N nodes. The original visibility graph (VG) algorithm implies that two arbitrary data values (t i , y i ) and (t j , y j ) have visibility, and consequently are connected nodes of the associated graph, if any other data (t k , y k ) placed between them (t i < t k < t j ) satisfies (Lacasa and Toral 2010): The schematic of the above geometric criterion and its associated visibility graph are shown in Fig. 1a. It is clear that two arbitrary nodes i and j in the graph are connected if one can draw a (straight) line in the time series joining y i and y j that does not intersect any intermediate data height.
The horizontal visibility graph (HVG) algorithm is an alternative, simpler and faster geometric criterion (Luque et al. 2009). In the HVG algorithm, two nodes i and j in the graph are connected if the following geometrical criteria is fulfilled within the time series: Figure 1b illustrates the scheme of the horizontal visibility algorithm and the associated graph. According to the HVG geometric criterion, two data points (t i , y i ) and (t j , y j ) are connected if one can draw a horizontal line in the time (2) y i , y j > y k , for all k such that: i < k < j series linking the corresponding imaginary vertical lines without intersecting any intermediate data point.
After constructing the original VG or HVG, the degree of each node is determined. The degree of node i is the number of links which touch i. The corresponding degree sequences (DS) of the original VG and the HVG algorithms are shown in Fig. 1c as time series. From Fig. 1c, it is clear that HVG is always a sub-graph of its associated original VG. This does not affect qualitative features of the resultant graph; quantitatively speaking, the HVG has less statistics (Luque et al. 2009). Now, the synchronization of the time series y 1 (t) and y 2 (t) is determined through computation of similarity of the DSs of the corresponding visibility graphs as: where DS (.) represents the degree sequence of a time series, and cov [.] and σ (.) are its covariance and the standard deviation. The synchronization values range from 0 to 1, where S = 0 means the time series are not synchronized, and S = 1 means that time series are identical.
The VG-based algorithms are nonlinear maps of time series to graphs, and it has been shown that many signal structural properties (e.g., periodicity, fractality, etc.) are inherited in the resultant graph (Yang et al. 2009;Korn and Faure 2003). However, some signal information is inevitably lost in the mapping from the fact that the network structure is completely determined by the binary adjacency matrix (Yang et al. 2009;Korn and Faure 2003). For example, two periodic signals with the same period would have the same visibility graph, albeit being quantitatively different. Furthermore, the VG algorithm is restricted to univariate time series, and since it is sensitive to the nonlinear monotonic transformation of the original time series, it cannot be used to represent invariants of the underlying system (Czechowski et al. 2016).
It is worth mentioning that the presented VG-based algorithms are computationally efficient to transform small-scale time series to graphs; however, it may take too much time to deal with very large time series. To transform a time series of size n, it is necessary check all the n(n-1)/2 pairs of signal nodes whether each pair of two nodes can see each other based on the defined geometrical criteria. For example, in Fig. 2a, to check the connection between signal nodes t 3 and t 7 , we need to know the maximum slope of the lines between node t 3 and nodes t 4 , t 5 and t 6 . After calculation for each pair of nodes, the maximum slope is timely updated. Hence, the time complexity of the connection judgment between each pair of nodes is O(1). Therefore, the total time complexity of the VG algorithm is O(n 2 ). In other words, the

State-space trajectory
Synchronization of chaos is often understood as a regime in which two coupled chaotic trajectories exhibit identical, but still chaotic oscillations. For example, the shapes of some well-known coupled chaotic systems seem regular, but the one-dimensional projection of their trajectories seems random. Many of the classical signal processing measures are not able solely to reveal such regularities behind the observed time series (Yang 2005;Packard et al. 1980). Therefore, a technique such as a state-space reconstruction is needed to detect this kind of behavior and analyze nonlinear time series. To this end, the time series of interest is embedded in a high-dimensional space to form a trajectory from which the properties of the original dynamic system can be deduced. The embedding theorem (or Taken's theorem Chan and Tong 2013) establishes that, when there is only a single sampled quantity from a dynamical system, it is possible to reconstruct a state-space that is equivalent to the original (but unknown) state-space composed of all the dynamical variables (Kantz and Schreiber 2004).
Two parameters time delay (T) and an embedding dimension (d) are defined to reconstruct a state-space. These parameters should be determined properly to avoid loss of information in the reconstructed state-space. It is suggested that to choose the same values of T and d for all signals to be able to compare the similarity of their states (Pereda et al. 2001). According to Taken's theorem (Chan and Tong 2013), the state vector of the i th sample of a given signal is defined as: Therefore, a multi-dimensional state-space is reconstructed from a scalar time series, as each trajectory contains [N − (d − 1)T] states, where N is the number of sampling points of each time series. (4) After mapping all samples of time series to the statespace, a window W n , with the width of 2(w 2 − w 1 ) is considered. Here, w 1 is the Theiler correction (Chan and Tong 2013) that is used to prevent information redundancy in the similarity computation, and w 2 is an integer number which determines the maximum temporal distance that a state can have from the reference state. Also, the state at the center of the window is considered as the reference state, Y c . Therefore, the width of window depicts the number of states which are restricted in the window, in other words, each window around reference state contains all states, Y m , as w 1 < |c − m| < w 2 . Now, by calculating the Euclidian distances of all states in the window with the reference state, the distance time series (DTS) of the window is constructed. After calculating the DTS, the original VG or the HVG algorithm is applied on the DTS to determine if the synchronization corresponds to the window W n . By shifting the window, the synchronization value is obtained for each window. Note that the window should be shifted with short enough steps to have a high temporal resolution which usually increases the computational cost. The overall synchronization of the system is obtained by averaging the computed synchronization values over all windows.
Using the Takens' theorem (Chan and Tong 2013), a signal sample is mapped nonlinearly to a state-space and the information of neighbor samples (i.e., past and next samples) are also used to create the state-space of that sample. The original VG and the HVG algorithms are also a nonlinear transform of signals to graphs and according to their geometric criteria, the information of at least nearest past and next samples is mapped into the graph. In other words, the original VG and the HVG algorithms inherently handle the high-dimensional chaos (Korn and Faure 2003), and we expect to observe no significant improvement on our analysis by making the trajectory of time series in state-space. However, in this work in addition to the original VG and the HVG algorithms, their combinations with the state-space are also examined as separate schemes, i.e. first, time series are mapped to the state-space and then the original VG (or HVG) algorithm is applied on the resultant DTS for each window. We call these schemes as the SS-VG and the SS-HVG, where SS refers to state-space mapping.

Results and discussions
We implemented the above VG-based similarity algorithms to find synchronization between chaotic (Sec. 3.1), noisy chaotic (Sec. 3.2) and stochastic (Sec. 3.3) time series. The performcance of the VG-based algorithms is compared against the commonly used synchronization measures: the crosscorrelation (CC), the coherence (Coh), the imaginary part of coherence (ImPC), the synchronization likelihood (SL), the phase coherence (PC) and the phase lag index (PLI). A brief review of these measures is provided in the Appendix.

Chaotic time series
Many output time series in physical, chemical and biological systems (e.g., EEG time seriess) have a chaotic motion in time. A chaotic motion means that the precise behavior of the system cannot be determined for a very long time in contrast to the periodic or quasi-periodic motion. Here, we use the coupled Rössler systems (Rössler 1976) to generate chaotic time series as (Smirnov and Andrzejak 2005): where the subscripts 1 and 2 denote the oscillator 1 (driver) with state vector of X = [x 1 , y 1 , z 1 ] and 2 (attractor) with state vector of Y = [x 2 , y 2 ,z 2 ], respectively. Parameters ω 1 = ω 0 − Δω and ω 2 + Δω are natural frequencies, where ω 0 = 1 is the normalized natural frequency and Δω is the natural frequency mismatch between two coupled systems. Two identical systems have the same parameters and therefore Δω = 0. It has been demonstrated that complete synchronization between two systems is not possible when there is a small but finite mismatch of the parameters of the systems (Osipov et al. 1997;Yanchuk et al. 2003).
The equations of systems are solved using the fourthorder Runge-Kutta method with a fixed step size of dt = 0.05. The initial conditions for all case studies are set as: x 1 (0) = 0.5, y 1 (0) = 1, z 1 (0) = 1.5, x 2 (0) = 2.5, y 2 (0) = 2 and z 2 (0) = 2.5. Each time series includes 10,000 data points, where the initial transients are removed by discarding the first 5000 data points. The coupling strength C is varied between 0 and 4, and the natural frequency mismatch of the oscillators is set as Δω = 0 and 0.05, which make the systems identical and non-identical, respectively. The different trajectories on the x 1 -x 2 plane obtained by changing the value of the coupling parameter (C) are illustrated in Fig. 2a-c when Δω = 0.05. In these figures, a clear tendency towards the identity of the two attractors can be observed, although complete synchronization will never be reached due to the parameter mismatch. In the Fig. 3 Predicted synchronization values as a function of coupling strength C. All compared measures are applied to time series x 1 and x 2 of two coupled Rössler systems case of complete synchronization, the x 1 -x 2 trajectory is presented as a straight line. Figure 3 shows the synchronization values computed by the original VG and the HVG measures as a function of coupling strength, compared to some common measures. The measures are applied on the x 1 and x 2 as a chaotic time series. The results for the identical Rössler system (i.e., Δω = 0) are presented in Fig. 3a. All the measures show an increase in the synchronization value as a function of C, and complete synchronization occurs for all measures roughly at C = 1.2. The results for a non-identical system with Δω = 0.05 are shown in Fig. 3b. For two nonidentical time series, the synchronization value should start with a value of approximately near zero and then increase gradually by increasing C. It can be seen that the predicted synchronization values by the HVG and SL measures are very close to 0 when C = 0, whereas other schemes start with larger values. All the measures ascend continuously regarding C, correctly representing a higher degree of synchronization when the coupling parameter is higher. All the measures, excluding the HVG and the PLI, approach to complete synchronization approximately at C = 2. The PLI scheme shows a constant value (roughly 0.55) for the synchronization even at very high coupling strength and does not show an ascending behavior. The HVG measure shows an ascending behavior over all values of the coupling strength and approaches slowly to complete synchronization. However, the synchronization values measured by the HVG never gets very close to full synchronization even at very high coupling strength (i.e., C = 25). This behavior is consistent with the concept of two non-identical systems which means that systems are not equal due to the mismatch parameter.
To further examine the accuracy and capability of the VG-based similarity measures, we define the synchronization error as (Yanchuk et al. 2003): Figure 4 shows x 1 (driver) and x 2 (respond) time series and corresponding errors time series of the identical coupled Rössler systems during an epoch of time when C = 0, 1 and 2, respectively.
It is clear from the figure that by increasing the coupling strength, the error decreases and two time series are matched and become synchronized. The amount of error after the transition region (i.e., time > 5000) around C = 1 is quite small. However, some error perturbations appeared for a while, which results in incomplete synchronization in this case. Consequently, these error perturbations affect the calculated synchronization values by the measures and as can be seen from Fig. 3a, some undulations appear around C = 1. As the coupling strength increases to values above 1.2, the amount of error becomes negligible, and measures should show complete synchronization.
Similarly, the driver and respond time series for the nonidentical coupled Rössler systems with Δω = 0.05 for C = 0, 2 and 5 along with corresponding errors are shown in Fig. 5. Due to the mismatch parameter, the error never becomes zero and as a results complete synchronization never is achieved. This figure shows that the amount of error around C = 0 is quite high and somehow is in the same order of magnitude with the time series. Therefore, we expect to have a weak coupling when C is close to zero. Therefore, we can say that the HVG and the (6) error = x 1 (t) − x 2 (t) Fig. 4 Time series depicting projective synchronization between drive state x 1 and response state x 2 (left) and the corresponding errors x 1 -x 2 of the identical coupled Rössler systems, C = {0, 1, 2} SL schemes predict the synchronization at weak coupling more accurate than the other measures (see Fig. 3b). Around C = 2, the value of the error is still significant, and we do not expect to have complete synchronization or values close to 1 for the synchronization value. Therefore, we can conclude that the PC, the CC, and the Coh measures overestimate the synchronization measures (see Fig. 3b). The behavior of the SL, the HVG, and the original VG measures are ascending as a function of coupling strength. However, due to the mismatch parameter, some errors always exist and two time series never match exactly. Figure 6 shows absolute value of the maximum errors for identical and non-identical coupled Rössler systems against coupling strength in a log-log plane. It can be seen that the slope of line for the identical case changes sharply at C = 1 and the error decreases dramatically. These changes in slope clearly explain the reason of complete synchronization around C = 1 predicted by various measures in Fig. 3a.
On the other hand, the slope of the line for the non-identical case is almost constant over all coupling strengths. It means that the error decreases monotonically as a function of C and becomes small only with very high C. Therefore, an accurate synchronization measure should also show a similar behavior for the non-identical systems, i.e., increasing gradually against increasing coupling strength. Therefore, we can say that the HVG and original VG measures, as well as the SL measure, show a good performance on finding the synchronization value for chaotic time series. However, the HVG measure shows a uniform increase after C = 1, whereas the VG and the SL schemes show a big jump in synchronization value and approach to complete synchronization when C is around 1 and then, show a uniform increase with a very smooth slope (see Fig. 3b). Our analysis shows that the maximum amount of error when C = 5 is about 10% of the magnitude of the time series, which means that the error is still significant. Therefore, we do not expect to have synchronization value close to 1 similar to predicted values by the original VG and the SL. Therefore, the original VG and the SL measures are quantitatively overestimating the synchronization compared to the HVG measure. Furthermore, due to the presence of the mismatch parameter, two time series never become synchronized completely. This concept is more consistent with the results of the HVG measure, which shows uniform increasing behavior over C with the almost constant increasing slope.
Note that the SS-VG and the SS-HVG overestimate the synchronization values for both the identical and non-identical cases. These schemes show the steepest increases (results are not shown here), whereas they almost predict complete synchronization for C < 0.5 for both the identical and non-identical Note that the computational parameters in this study for the state-space algorithm are time delay T = 2 and dimension d = 10. The small value of time delay helps to capture the shortest changes (e.g., high-frequency component) in the time series and increases the accuracy of the mapping. Also, since many real-world time series (e.g., the EEG or multisensor radar signals) have high intrinsic dimensionality, we selected d = 10 which would be a suitable value for our study.

Noisy chaotic time series
Output time series of many chaotic real-world systems are usually contaminated by noise. The signal-to-noise ratio (SNR) is a measure of signal strength relative to the level of noise in a system output and is usually measured in decibels (dB). To generate noisy chaotic time series, the Hénon map systems (Hénon 1976) are contaminated by additive white Gaussian noise (AWGN) with various SNRs (Diebold 1998). The Hénon map is defined as: (7) x 1,t+1 = 1 − 1.4x 2 1,t + by 1,t y 1,t+1 = x 1,t where the driver system includes states x 1 and y 1 , and the responder system contains states x 2 and y 2 . Here, C is the coupling parameter which defines the strength of the coupling between the two systems, and varies from 0 (indicating the subsystems are completely independent) to 1 (indicating the subsystems are completely synchronized). Also, t is discrete time or iteration index, and b is the bifurcation parameter, as for b = 0.3 both dynamic subsystems are identical, and for b ≠ 0.3 both dynamic subsystems are not identical (Sprott 2006). The different trajectories on the x 1 -x 2 plane for various values of the coupling parameter (C) for both noise-free and noisy identical cases are shown in Fig. 7a-c. A clear tendency towards the identity of the noise-free systems can be observed. However, for the noisy cases, the pattern of the response system is far from the driver system due to the presence of the noise and they never become matched even at C = 1. Figure 8 shows the performance of the VG-based synchronization measures on the noise-free identical and nonidentical Hénon map systems, compared to some common synchronization measures. The synchronization value should start with a value of approximately near zero when C = 0, and then increase gradually by increasing C. For both the identical and non-identical systems, all measures show an increase in synchronization value as the coupling strength increases. For the identical systems, a sharp transition in synchronization is observed almost for all measures when C becomes larger than 0.6, and complete synchronization occurs for all measures roughly at C > 0.7. One can see that all the VG-based measures exhibit good performances. These measures predict the synchronization value close to zero when C = 0, uniform increase up to C = 0.6 and finally complete synchronization for C > 0.7. Other common measures also show acceptable performance for the identical systems except for the Coh and the PC measures that overestimate the synchronization when there is a weak coupling between the systems (i.e., C < 0.5). The same story is observed for the non-identical systems. All the VG-based measures show an ascending behavior over all values of the coupling strength and approach gradually to 0.7 when C = 1. Note that for two non-identical systems, a complete synchronization is never achieved.
The results for the identical and non-identical noisy coupled Hénon systems are shown in Fig. 9 for SNR = 20 and 5 corresponds to low and high noise systems. The trend of all measures for the identical systems with a small amount of noise (i.e., SNR = 5, see Fig. 9a) are similar to the noise-free systems, i.e., by increasing the coupling strength, the measures should increase according to with the growing influence of the driver system on the response system. However, due to the presence of noise, some of the measures, including the VG-based measures, fail to predict complete synchronization at high coupling strengths (i.e., C > 0.7). It is clear that the Coh, PC, and CC measures predict complete synchronization at high coupling strengths which means that these measures are less sensitive to the noise at high coupling strengths. However, the Coh and PC measures overestimate the synchronization at low coupling strengths and predict large values for the synchronization even at C = 0. By increasing the amount of noise (i.e., SNR = 5. See Fig. 9b), all the measures fail to predict complete synchronization for the identical systems at high coupling strengths. The PC measure shows less sensitivity to the noise and predicts a higher value for the synchronization for larger coupling strengths. However, as we mentioned earlier, this measure fails to predict correct values at low coupling strengths. Among the VG-based techniques, the HVG measure shows less sensitivity to the noise and predicts a larger value for the synchronization at higher C compared to the other VG-based measures.
Since the non-identical systems are not similar, a complete synchronization is not achieved. Now, by adding some amount of noise to the non-identical coupled Hénon systems, the capability of the synchronization measures are examined. The results are shown in Fig. 9c, d for SNR = 20 and 5, respectively. Compared to the noise-free non-identical systems (see Fig. 8b), all the measures predict lower values for the synchronization due to the presence of the noise. The effect of noise become pronounced as the SNR decreases and consequently, a smaller value is predicted by all measures when SNR = 5 compared to SNR = 20.
To slightly go deeper on the effect of noise on the performance of the VG-based measures, we focus on the value of the synchronization when C = 1. For the noise-free nonidentical systems (see Fig. 8b), the synchronization value predicted by the HVG, the original VG, SS-HVG and SS-VG measures are equal to 0.8, 0.72, 0.8 and 0.75, respectively. At SNR = 20, the predicted values are 0.77, 0.63, 0.74 and 0.71 (see Fig. 10c) and at SNR = 5, the values are 0.42, 0.35, 0.3 and 0.25 (see Fig. 10d), respectively. These results indicate that the amount of noise when SNR = 20 leads to roughly 4, 13, 8 and 5% decrease in the predicted synchronization values by the HVG, original VG, SS-HVG and SS-VG measures compared to the noise-free systems, respectively. By decreasing the SNR to SNR = 5, the measures show more sensitivity to the noise and the differences increase to 47, 51, 62 and 67% for the HVG, the VG, the SS-HVG and the SS-VG measures, respectively. Therefore, we can say that the HVG measure shows less sensitivity to the noise compared to other VG-base measures. Also, we can conclude that the mapping the HVG and VG measures to the state-space leads to no improvement in the performance of the measures for the chaotic and noisy chaotic systems. Figure 10 shows absolute value of the maximum errors for (non-)identical noise-free/noisy coupled Hénon map systems against coupling strength in a semi-log plane. It is evident that by decreasing the SNR, the amount of noise and consequently the absolute error increases for both identical and non-identical systems. For the identical noise-free systems, a very sharp decrease in error is observed in C = 0.6, as the magnitude of error decreases five orders from C = 0.6 to C = 0.7. This sharp decrease in error means that two time series become quite synchronized for C > 0.6. By referring to Fig. 8a, we can see that all the measures show coordinated behavior and predict complete synchronization for C > 0.6. By adding noise to the identical systems, the smaller decrease in error is observed at C = 0.6, as the magnitude of error decreases only one order when SNR = 20. For the highly noisy system (i.e., SNR = 5) no changes in the order of error are observed. This large amount of noise in the systems contaminates the time series and the VG-based measures fail to distinguish accurately between real-time series and noise (see Fig. 9a, b). When coupling strength is high, the PC measure shows less sensitivity to noise and exhibits a good performance in detecting the synchronization of the noisy systems. A similar behavior can be observed for the non-identical systems. However, since the time series are not similar in non-identical systems, a complete synchronization is not achieved even for noise-free systems with large coupling strengths.
In the last part of this section, we study a quantity called graph entropy associated with the graph. There are several graph entropy calculation methods which incorporate random walks, degree distribution, and node centrality (Dehmer and Mowshowitz 2011). By calculating the entropy of the graphs corresponding to the noisy time series, the complexities of them compared to the noise-free case are measured. For this purpose, we use the so-called Shannon entropy formula (Shannon 2001) as: where p(k) is the probability of the degree k over the degree sequence. It is calculated by calculating the ratio of the number of nodes with degree k to the size of the degree sequence. The degree probabilities of the noise-free and noisy nonidentical Hénon map time series with C = 0 calculated by the HVG are shown in Fig. 11. It is clear that by contaminating the time series with noise, the fluctuations in the system and its DS become pronounced and consequently the probability of degrees with higher degrees increases. Needless to say, by increasing probability of higher degrees, the probability of lower degrees (e.g., k = 2) decreases.
We calculated the graph entropy for various values of coupling strength. As an example, the value of graph entropy, h, for the responding system of the non-identical Hénon map with C = 0 is equal to 2.47, 2.60 and 2.71 for noise-free, SNR = 20 and SNR = 5, respectively. We can see that by adding more noise to the system, the entropy of the system increases. Generally, the value of the entropy defines the uncertainty or randomness in the system and in a noisy system with an increase in noise, the entropy, and complexity of the system increases that lead to uncertainty about the information content of the time series. Therefore we can say that, the more noisy time series, the more complex graph and the more graph entropy. Note that this increase in entropy is observed for all values of coupling parameter, C, and also for the identical system as well.

Stochastic time series
To study the capability of the similarity measures for detecting real changes in synchronization of stochastic time series, we apply the measures on output time series of the so-called Kuramoto model (Kuramoto 2012). This model posits that the activity of a local stochastic system (e.g., neurons in the brain) can be sufficiently represented by its circular phase alone (Breakspear et al. 2010;Acebrón et al. 2005;Stam et al. 2007). The Kuramoto model describes the phase dynamics of a large network of N globally coupled limitcycle oscillators, θ i (t), having natural frequencies ω i (t), and whose phase dynamics is defined by the following differential equation: where K is the strength of the couplings between the oscillators. Actually, the i-th oscillator adjusts its phase velocity according to input from other oscillators through the coupling strength K. Usually, the global degree of synchrony of the system is described by an order parameter r(t) as: The averaged value of r(t) over time is abbreviated by r. The averaged order parameter r captures the degree of phase coherence in the system (Stam et al. 2007) as it approaches zero when the phases are uniformly distributed (have large Fig. 11 The node degree distribution of noise-free and noisy time series of the responding system of the Hénon map with C = 0 circular variance) and approaches one when the phases of all oscillators become perfectly locked, describes zero phase lag synchronization. The Kuramoto model shows a phase transition from a desynchronized to a partially synchronized state at a critical value K = K c . When K < K c , the system is not synchronized, and r = 0 (in the limit of an infinite number of oscillators, i.e., N → ∞ ). When K > K c , a single cluster of synchronized oscillators emerges, which grows for increasing K. For K > K c , the order parameter r is given by: Theoretically, an infinite number of oscillators is necessary for the analytical results to hold. However, it has been shown that with only 64 oscillators the model can explain various empirical results quite well (Stam et al. 2007;Kiss et al. 2002). Therefore, we solve the Kuramoto model for N = 64 in our work for various values of coupling strength K ranging from 0 to 8. Figure 12a-d show nil, weak, intermediate and full phase locking behavior of the selected number of oscillators for K = 0, 1, 3 and 8, respectively, in the polar form on a unit circle. The impact of increasing K on the model is increasing the phase synchrony amongst the oscillators. For K = 0 and 1, the oscillators disperse, whereas for K = 8 they are fully synchronized. For K = 3 we see that a large cluster of synchronous oscillators is apparent. However, some other oscillators, whose natural frequencies are at the tails of the distribution, are not locked to this cluster.
The Kuramoto model (Eq. 9) is solved numerically for each oscillator with a time step of 0.05 s. The simulation is run for 10,000 time steps and then the time series corresponding to oscillator i at time t is defined as x i = Sinθ i . It is clear that in the nil and weak cases (i.e., K = 0 and 1), the dynamic of time series are independent whereas for strong coupling (i.e., K = 8), the time series become relatively fully coherent. Thus, as the coupling strength K increases, the interaction functions to overcome the dispersion of natural frequencies ω i resulting in a transition from incoherence, to partial and then full synchronization.
For each value of K, five trials are done and the resulting time series of N = 64 oscillators are subjected to synchronization analysis. The final value of the synchronization at each K is computed by averaging the estimated synchronization over all possible pairs of the simulated oscillators (time series).
The values of the synchronization for the Kuramoto model calculated by the VG and the HVG similarity measures are shown in Fig. 13.
The results show the averaged values over all pairs of the 64 time series as a function of coupling strength. The results are compared with the averaged order parameter r (Eq. 10) as well as the results of the PLI, PC, and ImPC measures. It can be seen that r remains close to 0 until K reaches a critical value K c ≈ 1.5, above which r rapidly increases towards its asymptotic value of 1. A similar treatment is observed for the VG and the HVG similarity measures. However, compared to the r, the sharp increase happens in smaller K when using the original VG. The behavior of the HVG as well as the PLI and the PC measures are in close agreement with the results for r. The ImPC generally underestimated the true synchronization in the model. Its results starts to increase for K > K c , but never reached values much higher than 0.2 even for very high coupling strength K. Note that the SS-VG and the SS-HVG results (not shown here) show behavior similar to the VG and the HVG results and similar to the chaotic analysis, but they lead to no significant improvement in the results.
It is worth mentioning that the value of synchronization should be theoretically equal to zero for K < K c . The nonzero values below K c merely reflect fluctuations in the simulation due to a finite number of oscillators, N.

Application to real data
Our results showed that the HVG scheme is a capable and accurate method for finding synchronization between time series. For further study, we apply the HVG on real data.
It has been demonstrated that the degree distribution on the visibility graph associated with a random time series is presented theoretically by P(k) = (1/3)(2/3) (k−2) , where k is the degree of a node (Luque et al. 2009). Now, we apply the HVG on the data from the US Energy Information Administration on the crude oil future contract 1 (Dollars per Barrel) from Jan. 1986-Dec. 2000 (see Fig. 14a). A finite size time series is generated by extracting 200 sample data points of each year. Figure 14b shows a comparison between the results of HVG and the theory of data points extracted from the crude oil price. A fair agreement is observed for the degree distribution for small degrees. However, the method deviates from the theory for larger node degrees which are caused solely by finite size effects.
As a further study, the presented HVG-based synchronization method is applied on the brain EEG data to construct synchronization (or correlation) matrix of the brain. Electroencephalography or EEG is a well-known technique for measuring the brain activity. EEG signals have been considered to result either from random processes or to be generated by non-linear dynamic systems exhibiting chaotic behavior (Pijn et al. 1991). Furthermore, EEG time series can be contaminated with noise caused by user's head movements. We use the dataset that is described in (Hoffmann et al. 2008) for a healthy and a disabled subject (male, 51 years old) with multiple sclerosis (MS) disorder. The data is publicly available at http://mmspg .epfl.ch/BCI datas ets. The EEG was recorded at 2048 Hz sampling rate from 16 electrodes placed at the standard positions of the 10-20 international system (see Fig. 15). The indices of For the data described above, we build the adjacency matrices for the healthy and patient subjects. For this purpose, one needs to calculate the correlations among all pairs of EEG electrodes and deduce the respective adjacency matrix, called synchronization matrix. Therefore, the values of synchronization matrix describe the amount of synchronization between recording nodes. Hence, the presented HVG-based synchronization measure results in the calculation of a synchronization matrix which can be regarded as the adjacency matrix of a weighted graph (or network).
The synchronization matrix for the healthy subject is shown in Fig. 16a. A healthy brain network can be described as an intermediate between three extremes: a locally connected, highly ordered network; a random network; and a scale-free network (Stam 2014). Order is reflected in the high clustering of regular brain networks, whereas randomness (low order) is reflected in short path lengths. The scale-free component, characterized by a high diversity in node degree and high hierarchy, is indicated by the presence of highly connected hub nodes. The composite of these attributes in normal brains results in a hierarchial, modular network. It is clear that using the HVG-based method, a consistent network with the above definition for a healthy brain is achieved.
The synchronization matrix of the subject with MS disorder is shown in Fig. 16b. It has been shown that the functional connectivity decreases in the brain network of MS patients (Rocca et al. 2018). To better understand the effect of MS on the synchronization between brain channels, the synchronization difference between the healthy and patient subject is shown in Fig. 16c. Significant losses of synchronization are observed between Pz − P3, Pz − CP1, Pz − CP2, Pz − FC2, O1 − P3, C4 − P3, CP2 − Fz and CP2 − Cz. These results depicts that a significant MS lesion has been formed on the top of the brain. MS lesions cause segmental demyelination, axon energetic failure and axonal transection, which have a profound impact on nerve transmission and loss of synchrony by blocking nerve conduction.

Conclusion
The synchronization measures based on the visibility graph were used as a new way to detect synchronization between time series. By carrying out the VG algorithm, the data are transformed from the time domain to the graph (or network) domain, which allows for the dynamic of the time series to be studied via the topology of the graph network. In this work, the capabilities of the VGbased measures were analyzed by applying them to chaotic Rössler, noisy Hénon, and stochastic Kuramoto systems. For each case, the results were compared with some commonly used synchronization measures. Our study showed that the HVG and the VG similarity measures had a good performance in detecting the synchronization between chaotic and stochastic time series. However, compared to the original VG measure, the HVG similarity measure is more capable of detecting synchronization. Mapping the time series to the state-space and then applying the VG-based measures led to no improvement or no significant changes in the results. It means that the VG algorithm can capture the dynamic of high-dimensional time series quite well. On the other hand, the VG-based measures showed sensitivity to noise. They have reasonable performance on detecting synchronization with a low amount of noise; however, they did not have good capability in the high noisy systems. To summarize, the HVG similarity measure is reliable and robust in detecting synchronization in chaotic and stochastic time series. It also outperforms other common measures in detecting synchronization in noisy systems.
Our results are useful for choosing an adequate synchronization measure for different time series mining tasks including predictive modeling for diagnostics and pattern mining and clustering for knowledge discovery in different application domains including but not limited to brain informatics and medical informatics. Our further work includes using this results for detection of alcoholism, epilepsy and other disorders from EEG and functional brain network features extraction Pechenizkiy 2016, 2017).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.