Ordinal Synchronization: Using ordinal patterns to capture interdependencies between time series

To compute the (OS) between two time series $X$ and $Y$, we first extract their $D$-dimensional ordinal patterns [19]. In this way, we choose a length D and divide both time series of length M into $L=M/D$ equal segments. Next, we obtain the order of the values included inside each segment, also called the ordinal patterns: l

where $V_{t}$ and $W_{t}$ are the ordinal vectors inside the segment given by $\{t,t+1,...,t+D-1\}$, elements refer to the ordinal position of the values in $X_{t}$ and $Y_{t}$, respectively. Note that the elements in $V_{t}$ and $W_{t}$ are natural numbers ranging from $0$ to $D-1$. The higher the value in the time series, the higher the corresponding element in the ordinal vector. Following the example depicted in Fig. 1, where $D=4$, we obtain:

Then, we take the euclidean norm of each ordinal vector.

and we call $V^{N}_{t}=V_{t}/||V_{t}||$ and $W^{N}_{t}=W_{t}/||W_{t}||$ the normalized vectors. Note that this step only depends on the length D.

Now, we define the raw value of the instantaneous ordinal synchronization at time $t$ ($IOS^{raw}_{t}$) as the dot product between both ordinal ordinal vectors $IOS^{raw}_{t}=\sum_{i=1}^{D}v_{i,t}\textperiodcentered w_{i,t}$ (i.e., $V^{N}_{t}\cdot W^{N}_{t}$). For a more intuitive interpretation, we linearly rescale the value of $IOS^{raw}_{t}$ to be bounded between $-1$ and $1$:

where $min$ is the minimum possible value of the scalar product between two ordinal vectors. Note that, since the elements of the ordinal vectors are always positive and have only one component equal to zero, the lowest possible scalar product between $V_{t}$ and $W_{t}$ is obtained when the order of the elements of vector $V_{t}$ is inverted in $W_{t}$. In our example:

In general, for any vector of length D:

Following the normalization in 6), we ensure that two ordinal vectors that follow opposite evolutions will unambiguously lead to a value of $IOS_{t}=-1$, and two vectors whose elements have the same order will have an $IOS_{t}=1$. Being $L=M/D$ the total number of ordinal vectors in time series of $M$ points, the final value of the ordinal synchronization $OS\{X,Y\}$ for a given pair of time series $X$ and $Y$ is obtained averaging the instantaneous values of $IOS_{t}$ along the whole time series:

Since we consider the $IOS_{t}$ of consecutive (i.e., non-overlapping) time windows, the value of $t$ in Eq. 9 is given by the expression $t=1+iD$, with $i$ being a natural number bounded by $0\leq i\leq L-1$. Note that it is also possible to define a sliding $OS$ just by increasing $t$ in one unit for every $IOS_{t}$ instead of considering consecutive windows.

We analyzed the transition to the synchronized regime of two coupled Lorenz oscillators [20]. We implemented an electronic version of the Lorenz system, whose equations are detailed in Appendix B. Two Lorenz circuits are coupled unidirectionally in a master-slave configuration (see Fig. 2) with a coupling strength $\kappa$ that can be modified. Our experiments include two conditions: in the first one, $\kappa$ is modified in the absence of external noise; in the second one, $\kappa$ varies in presence of Gaussian noise with band selection. The (AI0-AI3) input ports of a data acquisition (DAQ) card are used for sampling the $x$ and $z$ variables of each circuit, while the output ports AO0 and AO1 generate two different noise signals ($\xi_{1}$, $\xi_{2}$) that perturb the dynamics of the Lorenz circuits through variable $x$ of each circuit. In this way, an external source of noise can be introduced to check the robustness of the experiments. The circuit responsible of the coupling strength $\kappa$ is controlled by a digital potentiometer XDCP, which is adjusted by digital pulses from ports P00 and P01. Noisy signals were designed in LabVIEW, using a Gaussian White Noise library [21] that generates two different Gaussian-distributed pseudorandom sequences bounded between [-1 1]. All the experimental process is controlled by a virtual interface in LabVIEW 2016 (PC).

The experiment works in the following way: First, $\kappa$ is set to zero and digital pulses (P00 and P01) are sent to the digital potentiometer until the highest value of $\kappa$ is reached. Second, variables $x$ and $z$ of the circuits are acquired by the analog ports (AI0-AI3) in order to compute the synchronization metrics. Initially, we have obtained all results for $\xi_{1}=\xi_{2}=0$, i.e., in the absence of external noise, and then, after a moderate amount of noise is introduced, all synchronization metrics are calculated again (See Appendix C). Every signal, with or without noise, has a length of 30000 samples.

We have checked the performance of the $OS$ in the context of neuroscientific datasets. Specifically, we quantified the level of synchronization between pairs of channels of MEG recordings. Data sets have been obtained from the Human Connectome Project (for details, see [22] and https://www.humanconnectome.org). The experimental data sets consist of $30$ MEG recordings of an individual during resting state for a period of approximately $2$ minutes each. During the scan, the subject were supine and maintained fixation on a projected red crosshair on a dark background. Brain activity was scanned with 241 magnetometers on a whole head MAGNES $3600$ (4D Neuroimaging, San Diego, CA, USA) system housed in a magnetically shielded room. The root-mean-squared noise of the magnetometers is about $5$ fT/sqrt ($Hz$) on average in the white-noise range (above 2 $Hz$). Data was recorded at sampling rate of $f_{s}\approx 508.63$ $Hz$. Five current coils attached to the subject, in combination with structural-imaging data and head-surface tracings, were used to localize the brain in geometric relation to the magnetometers and to monitor and partially correct for head movement during the MEG acquisition. Artifacts, bad channels, and bad segments were identified and removed from the MEG recordings, which were processed with a pipeline based on independent component analysis to identify and clean environmental and subject’s artifacts [22].

In order to assess the validity of $OS$, we have explored its performance for different values of $D$, from 3 to the full length of the time series under evaluation. Since it is the first time $OS$ is used, we have compared it to classical measures of correlation, namely Pearson correlation coefficient ($r$), spectral coherence ($SC$), phase locking value ($PLV$) and mutual information ($MI$). We have used two kinds of data sets to validate $OS$, on the one hand, experimental time series from nonlinear electronic circuits, and on the other hand, MEG recordings.

First, we take advantage of the ability of controlling the coupling strength between electronic circuits and investigate how $OS$ changes as two dynamical systems smoothly vary their level of synchronization from being unsynchronized to completely synchronized. Specifically, two electronic Lorenz systems are unidirectionally coupled with a parameter $\kappa$ controlling their coupling strength (see Appendix B for details). Initially, we do not perturb the oscillators with external noise (see Appendix C for the case of including external noisy signals). However, we can not avoid the intrinsic noise of the electronic circuits together with the tolerance of the electronic components (between 5 % and 10 %). Figure 3 shows how the value of $OS$ changes as the coupling strength $\kappa$ is increased from zero. Since the value of $OS$ depends on the length of the ordinal vectors, we show the results for three different values: $D=3$ (Fig. 3A), $D=500$ (Fig. 3B) and $D=1000$ (Fig. 3C). Note that, by increasing the length of the vectors, we are obtaining the amount of synchronization at different time scales. Together with $OS$, we plot the values of the rest of synchronization metrics in (A), (B) and (C), which remain unaltered in the three plots (since they do not depend on $D$).

In all cases, we observe that $OS$ increases for low to moderate values of $\kappa$ and remains at a high value once a certain threshold is reached. This behaviour is similar to the rest of the synchronization metrics. However, both $MI$ and $SC$ seem to saturate at values of $\kappa$ higher than $r$, $PLV$ and $OS$, which seem to reach a plateau around $\kappa=40$. Figure 3D shows the comparison of $OS$ for the three different values of $D$. Here, we can also observe how at $D=3$, $OS$ has a different qualitative behaviour from $D=500$ and $D=1000$, since it stays around 0.9 and does not reach 1 as in the windows of longer lengths. The reason is the existence of intrinsic noise of the electronic circuits, that affects much more the alignment of the ordinal vectors of shorter lengths than those with higher dimensions.

Figure 4 shows the average correlation ($\rho$) between each D-dependent $OS$ and the rest of synchronization metrics with zero noise. Note that correlations are higher than $0.92$ in all cases, although it seems to be certain vector lengths that maximize these correlations. Also note that correlations with $PLV$ and $r$ are the highest and, in all cases, very close to 1. At the same time, $MI$ and $SC$ show lower correlations that, in turn, seem to be more dependent on the value of the vector length $D$ (see Fig. 4A).

We can investigate how $OS$ is related to the rest of the synchronization metrics in more detail by setting the length of the ordinal vectors to a given value ($3$, $500$ or $1000$ in this case) and observe the influence of the level of synchronization (Fig. 4B, C and D). For any of the three selected lengths, $OS$ shows a linear relation with $PLV$ and $r$, especially at values of $OS$ higher than $0.5$. However, the relation with $SC$ and $MI$ seems to be nonlinear in all cases. Interestingly, for low levels of synchronization, $OS$ increases much faster than these two latter metrics. While $SC$ saturates around $0.8$, $MI$ finally increases faster than $OS$ only for high values of synchronization, eventually reaching the value of $OS$ around 1. Also note how, in the case of $D=3$ (Fig. 4B), the intrinsic noise of the electronic circuits prevents $OS$ to reach the value of one. This behaviour can be observed even clearer in the case of adding more noise into the system, as shown in Appendix C.

The second application is the evaluation of the level of synchronization between the 241 sensors measuring the activity of an individual during resting state. Concretely, we have 30 recordings of 2 minutes each. In this case, we can not control the amount of coupling between sensors but, alternatively, we have a diversity of levels of synchronizations between all possible pairs of sensors. Figure 5 shows how the correlations between $OS$ and the rest of the metrics change depending on $D$. As we can observe, correlations are high in all cases except for SC, but this one saturates around the same $D$ as the other synchronizations does.

As in the case of the electronic Lorenz oscillators, tuning the value of $D$ allows to obtain values of $OS$ closer, or more correlated, to other metrics. In fact, two different regions are clearly observed: (i) for values of $D\leq 20$ the correlation of $OS$ with PLV, r, and MI increases with $D$, while (ii) for $D>20$ correlation saturates around the highest value, being $r$ the metric with the highest correlation. Interestingly, the behaviour of the SC goes in the opposite direction, decreasing for higher values of $D$.

In order to gain insights about how the behaviour of $OS$ depends on the level of synchronization and the length $D$, we plot three different cases in Figure 6. In (A), we show the time series of two highly-correlated sensors, with their corresponding $OS$ value depending on $D$ (Fig. 6D). Plot (B) and (C) show the cases of two uncorrelated and negatively correlated sensors, respectively, with their values of $OS$ (Fig. 6E and F). Note that for the positive (negative) case, correlations tend to stabilize as $D$ grows, indicating the existence of a certain temporal scale at which synchronization is increased (reduced). Also note that, when time series are not correlated, this pattern is not that clear, and $OS$ values remain low for any value of $D$.

We had analyzed the relation of $OS$ with the rest of the metrics according to the level of synchronization. Figure 7 shows a panel of plots capturing the correlations between $OS$ and all other synchronization metrics for the MEG signals. Left plots show the case of $D=3$, middle plots show $D=500$ and right plots show $D=1000$. Different conclusions can be drawn depending on the synchronization metric $OS$ is compared to. In the case of $MI$ (first row), the existence of a nonlinear correlation between both metrics arises. However, this correlation decreases with the length of the ordinal vectors, becoming rather noisy for $D=3$. This behaviour is induced by the intrinsic noise of the MEG signals that, as in the case of electronic circuits, affects the value of $OS$ when short lengths of the ordinal vectors are considered. Also note that $MI$ is not able to distinguish positive from negative correlations between time series, a fact that makes $OS$ an interesting metric when both kind of synchronizations are expected. In our case, for example, despite the highest values of $OS$ are close to $1$, the lowest ones arrive to $-0.35$, indicating the existence of anti-correlated dynamics between certain pairs of sensors. A similar behaviour is reported in the case of the comparison with $PLV$ (second row). Again, a nonlinear relation exists between both metrics, which is rather noisy at low values of the ordinal vector lengths ($D=3$). $PLV$ has also the same limitations as $MI$, since it does not differentiate between positive and negative correlations. Interestingly, the relation with $SC$ is different from the two previous metrics (third row). Despite a nonlinear correlation between $OS$ and $SC$ seems to be present in the plots, this correlation is deteriorated with the increase of $D$.

Finally, $OS$ shows a clear linear correlation with $r$ (bottom row), which, as in the case of $MI$ and $PLV$ becomes noisy for low values of $D$. Note that, the loss of correlation for low values of $D$ is indicating that, at short time scales, $OS$ is capturing a different pattern of synchronization than at large scales. This is an interesting feature of $OS$ which suggests that, when using it as a metric to evaluate synchronization between signals, it is appealing to carry out an analysis depending on the vector length in order to reveal the existence of different levels of synchronization at different time scales.

The Pearson’s correlation coefficient $r$ consists of a covariance scaled by variances, thus capturing linear relationships among variables. From the equations of the variance (of $X$ and $Y$) and covariance (of $XY$), we obtain Pearson Correlation Coefficient as:

Pearson’s correlation is a measure of linear dependence between any pair of variables and it has the advantage of not requiring the knowledge of how variables are distributed. However, it should be applied only when variables are linearly related to each other.

Coherence (magnitude squared coherence or coherence spectrum) measures the linear correlation among the two spectra[12]. To calculate the coherence spectrum, data must be in the frequency domain. In order to do so, time series are usually divided into S sections of equal size. The Fast Fourier Transform algorithm is then computed over the sections to get the estimate of each section’s spectrum (periodogram). Then, the spectra of the sections is averaged to get the estimation of the whole data’s spectrum (Welch’s method). Finally, Coherence is a normalization of this estimate by the individual autospectral density functions [12]:

where $Sp_{xy}$ is the Cross Power Spectral Density (CPSD) of both signals, $Sp_{xx}$ and $Sp_{yy}$ are the Power Spectral Density (PSD) of the segmented signals $X$ and $Y$ taken individually, and $\langle\cdot\rangle$ is the average over the S segments. In the case of the data sets obtained with the nonlinear electronic Lorenz systems, frequencies higher than $f_{cut}=7.5K$ Hz have been disregarded for the computation of $SC$, since the power spectra of the electronic circuits are completely flat above this frequency. One of the drawbacks of Coherence is that it doesn’t discern the effects of amplitude and phase in the relationships measured between two signals, which makes its interpretation unclear [23, 4].

Phase Locking Value was first introduced by Lachaux et al. [23] as a new method to measure synchrony among neural populations. It has, at least, two major advantages over the classical coherence measure: it doesn’t require data to be stationary, a condition that can rarely be validated; and has a relatively easy interpretation (in terms of phase coupling). However, the methods used to extract instantaneous phase, a step needed to calculate $PLV$ rely on stationarity, so indirectly $PLV$ can be affected by this condition [24]. To obtain the $PLV$, the signal has to be decomposed to it’s instantaneous phases and amplitudes. To achieve this, there are several methods, such as Morlet wavelet convolution or Hilbert transform [12, 24]. In this work we will utilize the latter. Finally, $PLV$ is obtained averaging over time $t$:

where $\theta(t,n)$ is the (instantaneous) phase difference $\phi_{x}-\phi_{y}$, the phases to be compared from the signals $X$ and $Y$. Comparisons are carried out pairwise (bivariate).

Mutual Information is a measure of shared information between any components of a system, between systems, or any other parameter whose value’s probability can be estimated. It is based on Shannon’s notion of entropy, which, in a general sense, tries to quantify the amount of information contained in a random variable by means of its estimated probability distribution. Mutual information measures the amount of information shared between two random variables by means of its joint distribution, or conversely, the amount of information we can obtain from one random variable observing another. This is analogue to measuring the dependence between two random variables [25]. Let X and Y be two random variables with $\{{x_{1},x_{2},...x_{n}}\}$ and $\{{y_{1},y_{2},...y_{n}}\}$, $n$ possible values with probabilities $p(x)$ and $p(y)$. The $MI$ of X relative to Y can be written as:

where $p({x}\cap{y})$ is the probability that $X$ has a value of $x$ while $Y$ has a value of $y$, $H(X)$ is the entropy of $X$ and $H({X}|{Y})$ is the conditional entropy of $X$ and $Y$. One of the major advantages of $MI$ is that it captures linear and non-linear relationships among variables. One disadvantage is that it does not explicitly tell the shape of that distribution [24]. To get the mutual information between two random variables, we first need to estimate their probability density distribution [25, 24, 26]. Equation 16 compares joint probabilities against marginal ones. When two values are independent, the product of their marginal probabilities should equal their joint probability. When not, we can state that there is a relationship among them (not necessarily linear), because the probability of finding those values together is greater than the probability of finding them by chance. Thus, somehow, those time series are coupled, although we don’t know the way it occurs.