Assessing coupling dynamics from an ensemble of time series

Finding interdependency relations between (possibly multivariate) time series provides valuable knowledge about the processes that generate the signals. Information theory sets a natural framework for non-parametric measures of several classes of statistical dependencies. However, a reliable estimation from information-theoretic functionals is hampered when the dependency to be assessed is brief or evolves in time. Here, we show that these limitations can be overcome when we have access to an ensemble of independent repetitions of the time series. In particular, we gear a data-efficient estimator of probability densities to make use of the full structure of trial-based measures. By doing so, we can obtain time-resolved estimates for a family of entropy combinations (including mutual information, transfer entropy, and their conditional counterparts) which are more accurate than the simple average of individual estimates over trials. We show with simulated and real data that the proposed approach allows to recover the time-resolved dynamics of the coupling between different subsystems.

A technical problem that arises in various fields of science is to detect interdependencies between simultaneously measured time series. Namely, the detection of interdependencies is often the first step for elucidating how the subsystems that underly the time series interact. For example, in neuroscience, ecology, or econometrics this approach has lead to the discovery of new neural codes [1], better models of population dynamics [2], and methods to assess the influence of an economic variable [3], respectively. Tools to unveil an interdependency include linear techniques, such as cross-correlation and coherency analysis [4], non-linear synchrony measures [5], and the evaluation of statistical dependencies via mutual information (MI) [6]. Although useful for assesing the strength of the interaction, these indices do not allow to identify directionality (i.e. cause-response relationships). The latter are arguably more important, if one aims to understand the functioning of a system at the mechanistic level. In this paper we propose a method to reliably estimate the temporal course of directed interactions, as revealed by several information-theoretic functionals.
While causality is a broad concept, Wiener formulated an operative definition leaning on the idea that the cause occurs before the effect and, therefore, knowledge of the cause helps forecasting the effect [7]. The widely used Granger causality is the mathematical formalization of Wiener's definition in terms of linear regressions [4]. Alternatively, when a model of the underlying dynamics and of the interaction is not available, a sound nonparametric approach can be stated in terms of information theory. A prototypical example is transfer entropy (TE), which quantifies, in terms of a Kullback-Leibler divergence, how much the present and past of one system condition (i.e., cause in a Wiener sense) the dynamics of another [8]. Nevertheless, the non-parametric or modelfreeness nature of a measure does not usually come for free. A practical pitfall, common to most informationtheoretic approaches, is that they require many observations to be reliably estimated. This requisite directly confronts with situations in which the dependency to be analyzed evolves in time or is subjected to fast transients. When the non-stationarity is only due to a slow change of a parameter, over-embedding techniques can partially solve the problem by capturing the slow dynamics of the parameter as an additional variable [9]. It is also habitual to de-trend the time series or divide them into small windows within which the signals can be considered as approximately stationary. However, the above-mentioned procedures become unpractical when the relevant interactions change in a fast time scale. This is the common situation in brain responses and other complex systems where an external stimuli elicit a rapid functional reorganization of information-processing pathways.
Fortunately, in several disciplines the experiments leading to the multivariate time series can be systematically repeated. Thus, a typical experimental paradigm might render an ensemble of presumably independent repetitions or trials per experimental condition. In this letter we show that this multi-trial nature can be exploited to produce time-resolved estimates for a family of information-theoretic measures that we call entropy combinations. This family includes well-known functionals such as MI, TE, and their conditional counterparts: partial mutual information (PMI) [10] and partial transfer entropy (PTE) [11,12]. We use simulations and experimental data to demonstrate that the proposed ensemble estimators of entropy combinations are much more accurate than simple averaging of individual trial estimates.
We consider three simultaneously measured time series generated from stochastic processes X, Y , and Z which can be approximated as stationary Markov processes [13] of finite order. The state space of X can then be reconstructed using the delay embedded vectors x(n) = (x(n), ..., x(n − d x + 1)) for n = 1, . . . , N , where n is a time index and d x is the corresponding Markov order. Similarly we could construct y(n) and z(n) for processes Y and Z, respectively. Let V = (V 1 , ..., V m ) denote a random m-dimensional vector. Then, an entropy combination is defined by: [1,m] where χ S is the characteristic function of a set S. It can be easily checked that MI, TE, PMI and PTE are all entropy combinations: . The latter denotes the joint probability of finding X at states x(n + 1), x(n), ..., x(n − d x + 1) during time instants n + 1, n, n − 1, ..., n − d x + 1. Notice that, due to stationarity, p(x(n + 1), x(n)) is invariant under variations of the time index n.
A straightforward approach to the estimation of entropy combinations would be to add separate estimates of each of the involved multi-dimensional entropies. Popular estimators of differential entropy include plug-in estimators and fixed and adaptive histogram or partition methods. However, other non-parametric techniques such as kernel and nearest-neighbor estimators have been shown to be extremely more data-efficient and accurate [14]. An asymptotically unbiased estimator based on nearest-neighbor statistics is by Kozachenko and Leonenko (KL) [15]. For N realizations x[1], x [2], ..., x[N ] of a d-dimensional random vector X, the KL estimator takes the form: where ψ is the digamma function, v d is the volume of the d-dimensional unit ball, and ǫ(i) is the distance from x[i] to its kth nearest neighbor in the set {x[j]} ∀j =i . The KL estimator is based on the assumption that the density of the distribution of random vectors is constant within an ǫ-ball. The bias of the final entropy estimate depends on the validity of this assumption, and thus, on the values of ǫ(n). Since the size of the ǫ-balls depends directly on the dimensionality of the random vector, the biases of estimates for the differential entropies in (1) will, in general, not cancel, leading to a poor estimator of the entropy combination. This problem can be partially overcome by noticing that (2) holds for any value of k so that we do not need to have a fixed k. Therefore, we can vary the value of k in each data point so that the radius of the corresponding ǫ-balls would be approximately the same for the joint and the marginal spaces. This idea was originally proposed in [16] for estimating mutual information, was used in [17] to estimate PMI, and we generalize it here to the following estimator of entropy combinations: where F (k) = ψ(k) − ψ(N ) and · · · n = 1 N N n=1 (· · · ) denotes averaging with respect to the time index. The term k i (n) accounts for the number of neighbors of the nth realization of the marginal vector V Li located at a distance strictly less than ǫ(n), where ǫ(n) denotes the radius of the ǫ-ball in the joint space. The point itself is included in this counting.
A fundamental limitation of estimator (3) is the assumption that the involved multidimensional distributions are stationary. However, this is hardly the case in many real applications and time-adaptation becomes crucial in order to obtain meaningful estimates. A trivial solution is to use the following time-varying estimator of entropy combinations: This naive time-adaptive estimator is not useful in practice due to its large variance, which stems from the fact that a single data point is used for producing the estimate at each time instant. However, let us consider the case of an ensemble of r ′ repeated measurements (trials) from the dynamics of V . Let us also denote by v (r) [n] r the measured dynamics for those trials (r = 1, 2, ...r ′ ). Similarly, we denote by {v (r) i [n]} r the measured dynamics for the marginal vector V Li . A straightforward approach for integrating the information from different trials is to average together estimates obtained from individual trials: is the estimate obtained from the rth trial. However, this approach makes a poor use of the available data and will typically produce useless estimates, as will be shown in the experimental section of this chapter. A more effective procedure takes into account the multi-trial nature of our data by searching for neighbors across ensemble members, rather than n n* v (2) [n] n*+σ n*-σ (1) [n] n n* n*+σ n*-σ from within each individual trial. This nearest ensemble neighbors [18] approach is illustrated in Fig. 1 and leads to the following ensemble estimator of entropy combinations: where the counts of marginal neighbors {k are computed using overlapping time-windows of size 2σ, as shown in Fig. 1. For rapidly changing connectivity patterns, small values of σ might be needed to track the coupling dynamics while larger values of σ will lead to lower estimator variance.
To demonstrate thatΘ en can be used to characterize dynamic coupling patterns we simulated three nonlinearly coupled autoregressive processes with a timevarying coupling factor: during 1500 time steps and repeated R=50 trials with new initial conditions. The terms η x , η y and η z represent normally distributed noise processes, which are mutually independent across trials and time instants. The coupling delays amount to τ yx = 10, τ zy = 15 while the dynamics of the coupling follows a sinusoidal variation: Before PTE estimation each time-series was timedelayed so that they had maximal mutual information with the destination of the flow. That is, before computing some T a←b|c (n), the time-series b and c were delayed so that they shared maximum information with the timeseries a. To assess the statistical significance of the PTE values (at each time-instant) we applied a permutation test with surrogate data generated by randomly shuffling trials [19]. Fig. 2 shows the time-varying PTEs obtained for these data with the ensemble estimator of entropy combinations given in Eq. 5. Indeed, the PTE analysis accurately describes the underlying interaction dynamics. In particular, it captures both the onset/offset and the oscillatory profile of the effective coupling across the three processes. On the other hand, the naive average estimator (5) did not reveal any significant flow of information between the three time-series (see supplementary material).
To evaluate the robustness and performance of the entropy combination estimator to real levels of noise and measurements variability, we also present a second example derived from experimental data on electronic circuits. The system consists of two nonlinear Mackey-Glass circuits unidirectionally coupled through their voltage variables. The master circuit is additionally subject to a feedback loop responsible for generating high dimensional chaotic dynamics. A time-varying effective coupling is then induced by periodically modulating the strength of the coupling between circuits as controlled by an external CPU. In this case, we applied transfer entropy between the voltage signals generated from the two circuits for 180 trials, each 1000 sampling times long. Figure 3 shows the TE estimates obtained with (5) versus the temporal lag introduced between the two voltage signals (intended to scan the unknown coupling delay). The results show that the TE estimates capture perfectly the dynamics of the effect exerted by the master circuit on the slave circuit. On the other hand, no significant coupling is detected in the reverse direction (see supplementary material). Both the period of the coupling dynamics (100 samples) and the coupling delay (20 samples) can be accurately recovered from Fig. 3.
In conclusion, we have introduced an ensemble estimator of entropy combinations that is able to detect time-varying information flow between dynamical systems, provided that an ensemble of repeated measurements is available for each system. The proposed approach allows to construct time-adaptive estimators of MI, PMI, TE and PTE, which are the most common information-theoretic measures for dynamical coupling analyses. Using simulations and real physical measurements from electronic circuits we showed that these new estimators can accurately describe multivariate coupling dynamics. It is important to mention that intrinsic to our approach is the assumption that the evolution of the interdependencies to be detected are to some degree "locked" to the trial onset. This is typically the case when some controlled external perturbation induce or evokes the interactions across the subsystems measured. The degree of locking determines the maximum temporal resolution achievable by the method (which is controlled via σ). Nevertheless, alignment techniques can help to reduce the possible jitter across trials and thus increase the resolution. The methods presented here are general but we anticipate that a potential application is the analysis of the mechanisms underlying the generation of event-related brain responses and the seasonal variations of geophysical variables. To promote dissemination we have publicly released a software library that includes efficient implementations of these and other informationtheoretic methods [20]. This work has been supported by the EU project GABA (FP6-2005-NEST-Path 043309) and by the Finnish Foundation for Technology Promotion.