Experiments

Simulated Data

We simulated five time series in two cases: linear Gaussianly distributed, and non linear, generated by an AR model, equations (1), (2), Baccala (2001). The following equations are for the linear Gaussian autoregressive model:

(1)   \begin{equation*} \begin{aligned} 	X_{1,n} &=  0.95 \sqrt{2} \, X_{1,n-1} - 0.9025\, X_{1,n-2} + w_{1,n} \\ 	X_{2,n} &= 0.5\, x_{1,n-2} + w_{2,n}  \\ 	X_{3,n} &= -0.4\, X_{1,n-3} + w_{3,n}  \\ 	X_{4,n} &= -0.5\, X_{1,n-2} + 0.25 \sqrt{2} \, X_{4,n-1} + 0.25 \sqrt{2} \, X_{5,n-1} + w_{4,n}  \\ 	X_{5,n} &= -0.25 \sqrt{2} \, X_{4,n-1} + 0.25 \sqrt{2} \, X_{5,n-1}  + w_{5,n} \end{aligned} \end{equation*}

where w_{1,n}, w_{2,n}, w_{3,n}, w_{4,n}, w_{5,n} are drawn from Gaussian noise with zero mean and unit variance. The following are the equations for the non linear model:

(2)   \begin{equation*} \begin{aligned} 	X_{1,n} &=  0.95 \sqrt{2} \, X_{1,n-1} - 0.9025\, X_{1,n-2} + z_{1,n} \\ 	X_{2,n} &= 0.5\, X_{1,n-2}^2 + z_{2,n}  \\ 	X_{3,n} &= -0.4\, X_{1,n-3} + z_{3,n}  \\ 	X_{4,n} &= -0.5\, X_{1,n-2}^2 + 0.25 \sqrt{2} \, X_{4,n-1} + 0.25 \sqrt{2} \, X_{5,n-1} + z_{4,n}  \\ 	X_{5,n} &= -0.25 \sqrt{2} \, X_{4,n-1} + 0.25 \sqrt{2} \, X_{5,n-1} + z_{5,n}   \end{aligned} \end{equation*}

where z_{1,n}, z_{2,n}, z_{3,n}, z_{4,n}, z_{5,n} are drawn from Gaussian noise with zero mean and unit variance. A schematic representation of the simulated couplings, valid for both systems, is reproduced in “Simulated system” figure shown below.

Then we considered a challenging situation in terms of number of interacting systems and lag of the interaction effects, considering the time series simulated with equations (1) and (2), which involve five systems and contain influences up to 3 points in the past. The experiments were run on 100 realizations of eqs. (1) and (2), of length equal to 512 points. We investigated the TE between each pair of variables conditioned to the other three. The setup of the experiment was the following: for all estimators, used either in the UE or in the NUE framework, the maximum lag for the candidates was set as 5, the number of surrogates was fixed to 100 and \alpha = 0.05. We set 6 quantization levels for BIN and 10 nearest neighbours for NN estimator.

In order to check whether the methods were able to detect the right information transfers, taking into account figure “Symulated system”, we expected that the estimators found a TE greater than zero with the highest significance at the following matrix elements: (1,2), (1,3), (1,4), (4,5), (5,4). Figures “Results on linear AR” and “Results on non-linear AR” report the analysis results obtained respectively for the linear system and the non linear system. Looking at Figure “Results on linear AR” one can notice that LIN UE has very good performances: this reflects the fact that this approach is, “by construction”, the best way to detect information flows. Its NUE version can detect the same links between the variables, though with a slightly higher number of false positives. The LIN estimator, therefore, is able to reveal the correct information flows for this simulation. On the contrary, the BIN UE suffers from the curse of dimensionality mentioned previously, in Entropy estimators section: evaluating the influences up to the first 5 past points for all the series implies that the uniform embedding is projecting data into a phase space of M \times 5 dimensions, where M is the number of time series, so the phase space will have 25 dimensions and the points are spread enough to loose relevant information about the transfer entropies governing the system. Therefore, as we can see, BIN UE can not detect any expected link. NN UE retrieves all the true links, but also detects a non-negligible number of false interactions. On the other hand, BIN and NN used in the NUE framework are able to recover all the correct links, with only a few false links. Moving to Figure “Results on non-linear AR” depicting TE analysis for the non-linear systems, one can notice that the LIN estimator cannot detect all the correct information flows and returns also some false links. Again, BIN UE cannot detect any link because of the curse of dimensionality, while BIN NUE has high specificity and sensitivity because it can retrieve all the correct information transfers with limited detection of detection of false positives. NN NUE can obtain almost the same performance as BIN NUE but its specificity is lower, especially along the direction X_2 \rightarrow X_4. NN UE this time is not able to detect all the correct information transfers (X_5 \rightarrow X_4 remains undetected) and reveals some false links (X_2 \rightarrow X_4, X_2 \rightarrow X_5).

To better clarify whether and how much the methods are able to distinguish between the true information transfer links and the false ones, in Figures “TE vs Significance on linear AR” and “TE vs Significance on non-linear AR” we plotted the average TE with respect to the number of significant realizations found by the methods. Each retrieved link is a point in this bidimensional space. The true links should be in the upper right corner of the plot corresponding to high TE and high number of significant realizations, and they should be apart from the false links, whose natural location would be around the origin of the plot (low TE and low number of significant links). Looking at figure “TE vs Significance on linear AR” one can notice that for all the methods, except BIN UE and partly NN UE, the two groups of links are well separated and the false links with an averaged TE greater than zero in figure “Results on linear AR” can be neglected. The opposite reasoning holds for BIN UE that is not able to distinguish between false and true links. For the non linear system (figure “TE vs Significance on linear AR”) only BIN NUE can separate well true positive from false positive links.

To understand how stable the performance of the methods is, in terms of sensitivity and the specificity, with respect to the length of the analyzed data set, we computed the analysis varying the series length from 128 to 1024 points. Figures “ROC curve for linear AR” and “ROC curve for non-linear AR” depict, respectively for the systems (1) and (2), the ROC curves obtained for all methods as a function of the series length. Evaluating the amount of TP (true positives), TN (true negatives), FP (false positives) and FN (false negatives) after grouping together all coupled directions (positives) and all uncoupled directions (negatives), we computed sensitivity as TP / (TP + FN) and specificity as TN / (TN + FP). In the case of the linear system (Figure “ROC curve for linear AR”), all methods except the BIN UE provide good performance, with the LIN estimator providing the best sensitivity and specificity. All methods provided robust results robust with respect to the series length, with only a limited degrade in the performance observed for 128 points. In the case of the non-linear system (Figure “ROC curve for non-linear AR”), the performance was optimal for BIN NUE and NN NUE (with a slightly lower specificity), while the methods implementing either the LIN estimator or the UE approach were considerably less sensitive.

  • Example
  • Example
  • Example
  • Example
  • Example
  • Example
  • Example

Application to Cardio-Pulmonary system

We considered cardiorespiratory time series measured from 15 young healthy subjects (25.7 \pm 2.7 years old) undergoing a standard head-up tilt testing protocol, Faes (2011). The acquired signals were the surface ECG, the finger arterial blood pressure, and the respiratory nasal flow, measured at 1 kHz sampling rate for 15 minutes in the
resting supine position, and 15 further minutes in the 60^\circ position after passive head-up tilting of the bed table. From these signals, the beat-to-beat variability series of heart period (RR interval), RR(n), systolic arterial pressure (SAP), Sap(n), and respiratory activity, Resp(n), were offline measured respectively as the temporal interval occurring between the n-th and the (n+1)-th R waves of the ECG, as the local maximum of the systolic arterial pressure signal inside the n-th heartbeat, and as the nasal flow taken at the onset of the n-th heartbeat. This measurement convention allows instantaneous effects from Sap(n) to RR(n), as well as from Resp(n) to Sap(n) and to RR(n), which were implemented using the relevant feature of the toolbox. The subsequent data analysis was performed on stationary windows of 300 beats taken in supine and upright body positions; inside these windows, the series were normalized to zero mean and unit variance, obtaining the dimensionless series resp(n), sap(n), rr(n).

The analysis of the information transfer for cardiovascular and cardiorespiratory time series was focused on the directions of interaction that are more studied from a physiological point of view: the link from SAP to RR which is related to the so-called cardiac baroreflex, and the links originating from Resp and directed either to RR or to SAP, related respectively to cardiopulmonary or vasculo-pulmonary regulation mechanisms, Faes (2011). Figure “Supine vs Upright” reports the distribution of the multivariate TE computed along these directions using all methods, with subjects studied in the supine and upright body positions. In the studied protocol, the transition from supine to upright is known to evoke an activation of the sympathetic nervous system and a concurrent deactivation of the parasympathetic nervous system, Cohen (2002). Accordingly, the two main physiological regulation mechanisms that are expected to be solicited by this transition are: (i) a substantial increase of the baroreflex regulation (direction sap \rightarrow rr), reflecting the necessity of the cardiovascular system to react with changes in the heart rate to the higher fluctuations in the arterial pressure induced by the sympathetic activation; and (ii) a substantial decrease of cardiopulmonary regulation (direction resp \rightarrow rr), reflecting the dampening of respiratory sinus arrhythmia consequent to the parasymphatetic deactivation, Nollo (2009). On the contrary, no known alterations of the vasculo-pulmonary regulation (direction resp \rightarrow sap) are expected when moving from supine to upright, Faes (2011). In our analysis all these trends are well reflected in terms of information transfer when the multivariate TE is estimated using the BIN NUE and NN NUE methods. Indeed we observe in Figure “Supine vs Upright” that BIN NUE and NN NUE reveal, moving from supine to upright, a substantial increase of the TE from Sap to RR, a substantial decrease of the TE from Resp to RR, and an unchanged TE from resp to SAP. These trends were also observed, though with less evident differences, computing the TE according to the LIN estimator. These results suggest the appropriateness of model free TE estimators based on NUE for detecting the information transfer in physiological time series. On the contrary, the BIN UE estimator shows different trends of difficult physiological interpretation, thus suggesting also in experimental data that the estimated TE may be unreliable due to the curse of dimensionality.

  • Example

EEG Data

The following experiment is performed on intracranial electroencephalography (EEG) measurements recorded from a patient with refractory epilepsy. The dataset consists of time series from 76 contacts. The first sixty-four of these contacts were placed on a 8×8 grid at the cortical level, while the other 12 were along two six-electrode strips that were implanted in deeper brain structures. Eight sets of measurements were taken on this patient, corresponding to eight different epileptic seizures. An epileptologist, examining the data for each seizure, identified two key periods relating to the seizure i.e., a pre-ictal period, just before the clinical onset, and an ictal one, corresponding to the seizure spread and to the clinical symptoms. Each epoch contained 10 seconds of data recorded at 400 Hz. The data are available here and described in Kramer (2008). In order to reduce overfitting, data were downsampled to 100 Hz.

In such high dimensional and redundant data, a non-uniform embedding approach is by construction the most appropriate to identify the patterns of information transfer specific to the onset and spread of the epileptic seizure. The aim of the experiment was to use the NUE approach in order to try to distinguish the system behaviour looking at the information transferred between the variables during the pre-ictal and ictal phases. The embedding size was set to eight. The regions corresponding to one of the depth strips (contacts 70 to 76) and the lower left corner of the grid (contacts 1-4, 9-11 and 17) were resected during anterior temporal lobectomy as they were identified by the epileptologists as the seizure onset zone. The Binning approach to NUE seems to be the one which best identifies these areas as those most influential at the start of the seizure and in the early phases of the spread, signature of a putative seizure onset zone. The Binning approach is more selective with respect to the target variables for each driver and less sensitive to conduction effects resulting in the diagonal patterns observed with the other methods and probably done to conduction effects on the grid.

  • Example

Conclusions

I presented three entropy estimators able to reveal the information transferred among variables represented by time series. I implemented the estimators in two different ways according to UE and NUE approaches, resulting in six methods, two of which are novel, BIN NUE and NN NUE. I compared all the methods validating them on simulated data first and then on real data. I checked whether and how the methods were affected by the number of variables and by the time lag at which the series influenced each other. From the results obtained we can conclude that the new methods introduced, not assuming any model to explain the data and exploiting the NUE strategy for component selection, can detect the correct information flows and are less affected by the number of involved processes and by their interaction lags. The NUE approaches are indeed prone to work in high dimensional spaces as well as in low dimensional spaces because of their ability to reduce the effective dimension of the phase space, choosing only the right variables at the specific time lag that are better able to explain the destination series. On the contrary, BIN UE and NN UE suffer from the curse of dimensionality when several time series and longer interaction delays are present. Finally, looking at LIN UE and LIN NUE performances we can conclude that, even though the equivalence between Granger causality and TE establishes a convenient joint framework for these two measures, there are some drawbacks in having a predefined model to explain the data when these are non linear. The better performances obtained by the new methods appear when looking at the ROC curves: BIN NUE and NN NUE have high sensitivity and specificity both for linear and non linear systems.

Bibliography

  1. Baccal'a, Luiz A, Sameshima, Koichi: Partial directed coherence: a new concept in neural structure determination. In: Biol Cybern, 84 (6), pp. 463–474, 2001.
  2. Marinazzo, Daniele, Pellicoro, Mario, Stramaglia, Sebastiano: Kernel method for nonlinear Granger causality. In: Phys Rev Lett, 100 (14), pp. 144103, 2008.
  3. Faes, Luca, Nollo, Giandomenico, Porta, Alberto: Information domain approach to the investigation of cardio-vascular, cardio-pulmonary and vasculo-pulmonary causal couplings. In: Front Physiol, 2 (80), pp. 80, 0000.
  4. Cohen, Michael A, Taylor, J Andrew: Short-term cardiovascular oscillations in man: measuring and modelling the physiologies. In: Physiol J, 542 (3), pp. 669–683, 2002.
  5. Nollo, Giandomenico, Faes, Luca, Antolini, Renzo, Porta, Alberto: Assessing causality in normal and impaired short-term cardiovascular regulation via nonlinear prediction methods. In: Philos Trans A Math Phys Eng Sci, 367 (1892), pp. 1423–1440, 2009.
  6. Kramer, Mark A, Kolaczyk, Eric D, Kirsch, Heidi E: Emergent network topology at seizure onset in humans. In: Epilepsy Res, 79 (2), pp. 173–186, 2008.