Multiscale Information Decomposition: Exact Computation for Multivariate Gaussian Processes

Abstract

Exploiting the theory of state space models, we derive the exact expressions of the information transfer, as well as redundant and synergistic transfer, for coupled Gaussian processes observed at multiple temporal scales. All of the terms, constituting the frameworks known as interaction information decomposition and partial information decomposition, can thus be analytically obtained for different time scales from the parameters of the VAR model that fits the processes. We report the application of the proposed methodology firstly to benchmark Gaussian systems, showing that this class of systems may generate patterns of information decomposition characterized by prevalently redundant or synergistic information transfer persisting across multiple time scales or even by the alternating prevalence of redundant and synergistic source interaction depending on the time scale. Then, we apply our method to an important topic in neuroscience, i.e., the detection of causal interactions in human epilepsy networks, for which we show the relevance of partial information decomposition to the detection of multiscale information transfer spreading from the seizure onset zone.

1. Introduction

The information-theoretic treatment of groups of correlated degrees of freedom can reveal their functional roles as memory structures or information processing units. A large body of recent work has shown how the general concept of “information processing” in a network of multiple interacting dynamical systems described by multivariate stochastic processes can be dissected into basic elements of computation defined within the so-called framework of information dynamics [1]. These elements essentially reflect the new information produced at each moment in time about a target system in the network [2], the information stored in the target system [3,4], the information transferred to it from the other connected systems [5,6] and the modification of the information flowing from multiple source systems to the target [7,8]. The measures of information dynamics have gained more and more importance in both theoretical and applicative studies in several fields of science [9,10,11,12,13,14,15,16,17,18]. While the information-theoretic approaches to the definition and quantification of new information, information storage and information transfer are well understood and widely accepted, the problem of defining, interpreting and using measures of information modification has not been fully addressed in the literature.

Information modification in a network is tightly related to the concepts of redundancy and synergy between source systems sharing information about a target system, which refer to the existence of common information about the target that can be retrieved when the sources are used separately (redundancy) or when they are used jointly (synergy) [19]. Classical multivariate entropy-based approaches refer to the interaction information decomposition (IID), which reflects information modification through the balance between redundant and synergetic interaction among different source systems influencing the target [20,21,22]. The IID framework has the drawback that it implicitly considers redundancy and synergy as mutually exclusive concepts, because it quantifies information modification with a single measure of interaction information [23] (also called co-information [24]) that takes positive or negative values depending on whether the net interaction between the sources is synergistic or redundant. This limitation has been overcome by the elegant mathematical framework introduced by Williams and Beer [25], who proposed the so-called partial information decomposition (PID) as a nonnegative decomposition of the information shared between a target and a set of sources into terms quantifying separately unique, redundant and synergistic contributions. However, the PID framework has the drawback that the terms composing the PID cannot be obtained unequivocally from classic measures of information theory (i.e., entropy and mutual information), but a new definition of either redundant, synergistic or unique information needs to be provided to implement the decomposition. Accordingly, much effort has focused on finding the most proper measures to define the components of the PID, with alternative proposals defining new measures of redundancy [25,26], synergy [27,28] or unique information [29]. The proliferation of different definitions is mainly due to the fact that there is no full consensus on which axioms should be stated to impose desirable properties for the PID measures. An additional problem which so far has seriously limited the practical implementation of these concepts is the difficulty in providing reliable estimates of the information measures appearing in the IID and PID decompositions. The naive estimation of probabilities by histogram-based methods followed by the use of plug-in estimators leads to serious bias problems [30,31]. While the use of binless density estimators [32] and the adoption of schemes for dimensionality reduction [33,34] have been shown to improve the reliability of estimates of information storage and transfer [35], the effectiveness of these approaches for the computation of measures of information modification has not been demonstrated yet. Interestingly, both the problems of defining appropriate PID measures and of reliably estimating these measures from data are much alleviated if one assumes that the observed variables have a joint Gaussian distribution. Indeed, in such a case, recent studies have proven the equivalence between most of the proposed redundancy measures to be used in the PID [36] and have provided closed form solutions to the issue of computing any measure of information dynamics from the parameters of the vector autoregressive (VAR) model that characterizes an observed multivariate Gaussian process [17,37,38].

The second fundamental question that is addressed in this study is relevant to the computation of information dynamics for stochastic processes displaying multiscale dynamical structures. It is indeed well known that many complex physical and biological systems exhibit peculiar oscillatory activities, which are deployed across multiple temporal scales [39,40,41]. The most common way to investigate such activities is to resample at different scales, typically through low pass filtering and downsampling [42,43], the originally measured realization of an observed process, so as to yield a set of rescaled time series, which are then analyzed employing different dynamical measures. This approach is well established and widely used for the multiscale entropy analysis of individual time series measured from scalar stochastic processes. However, its extension to the investigation of the multiscale structure of the information transfer among coupled processes is complicated by theoretical and practical issues [44,45]. Theoretically, the procedure of rescaling alters the causal interactions between lagged components of the processes in a way that is not fully understood and, if not properly performed, may alter the temporal relations between processes and thus induce spurious detection of information transfer. In practical analysis, filtering and downsampling are known to degrade severely the estimation of information dynamics and to impact consistently the detectability, accuracy and data demand [46,47].

In recent works, we have started tackling the above problems within the framework of linear VAR modeling of multivariate Gaussian processes, with the focus on the multiscale computation of information storage and information transfer [48,49]. In this study, we aim at extending these recent theoretical advances to the multiscale analysis of information modification in multivariate Gaussian systems performed through the IID and PID decomposition frameworks. To this end, we exploit the theory of state space (SS) models [50] and build on recent theoretical results [44,45] to show that exact values of interaction transfer, as well as redundant and synergistic transfer can be obtained for coupled Gaussian processes observed at different time scales starting from the parameters of the VAR model that fits the processes and from the scale factor. The theoretical derivations are first used in examples of benchmark Gaussian systems, reporting that these systems may generate patterns of information decomposition characterized by prevalently redundant or synergistic information transfer persisting across multiple time scales or even by alternating the prevalence of redundant and synergistic source interaction depending on the time scale. The high computational reliability of the SS approach is then exploited in the analysis of real data by the application to a topic of great interest in neuroscience, i.e., the detection of information transfer in epilepsy networks.

The proposed framework is implemented in the msID MATLAB^® toolbox, which is uploaded as Supplementary Material to this article and is freely available for download from www.lucafaes.net/msID.html and https://github.com/danielemarinazzo/multiscalePID.

2. Information Transfer Decomposition in Multivariate Processes

Let us consider a discrete-time, stationary vector stochastic process composed of M real-valued zero-mean scalar processes, ${\mathit{Y}}_n={[Y_{1,n}\cdots Y_{M,n}]}^T$, $-\infty <n<\infty$. In an information-theoretic framework, the information transfer between scalar sub-processes is quantified by the well-known transfer entropy (TE), which is a popular measure of the “information transfer” directed towards an assigned target process from one or more source processes. Specifically, the TE quantifies the amount of information that the past of the source provides about the present of the target over and above the information already provided by the past of the target itself [5]. Taking $Y_j$ as target and $Y_i$ as source, the TE is defined as:

$${\mathcal{T}}_{i\to j}=I(Y_{j,n};Y_{i,n}^{-}|Y_{j,n}^{-})$$

(1)

where $Y_{i,n}^{-}=[Y_{i,n-1},Y_{i,n-2}\cdots ]$ and $Y_{j,n}^{-}=[Y_{j,n-1},Y_{j,n-2}\cdots ]$ represent the past of the source and target processes and $I(·;·|·)$ denotes conditional mutual information (MI). In the presence of two sources $Y_i$ and $Y_k$ and a target $Y_j$, the information transferred toward $Y_j$ from the sources $Y_i$ and $Y_k$ taken together is quantified by the joint TE:

$${\mathcal{T}}_{ik\to j}=I(Y_{j,n};Y_{i,n}^{-},Y_{k,n}^{-}|Y_{j,n}^{-}).$$

(2)

Under the premise that the information jointly transferred to the target by the two sources is different than the sum of the amounts of information transferred individually, in the following, we present two possible strategies to decompose the joint TE (2) into amounts eliciting the individual TEs, as well as redundant and/or synergistic TE terms.

2.1. Interaction Information Decomposition

The first strategy, which we denote as interaction information decomposition (IID), decomposes the joint TE (2) as:

$${\mathcal{T}}_{ik\to j}={\mathcal{T}}_{i\to j}+{\mathcal{T}}_{k\to j}+{\mathcal{I}}_{ik\to j},$$

(3)

where $I_{ik\to j}$ is denoted as interaction transfer entropy (ITE) because it is equivalent to the interaction information [23] computed between the present of the target and the past of the two sources, conditioned to the past of the target:

$${\mathcal{I}}_{ik\to j}=I(Y_{j,n};Y_{i,n}^{-};Y_{k,n}^{-}|Y_{j,n}^{-}).$$

(4)

The interaction TE quantifies the modification of the information transferred from the source processes $Y_i$ and $Y_k$ to the target $Y_j$, being positive when $Y_i$ and $Y_k$ cooperate in a synergistic way and negative when they act redundantly. This interpretation is evident from the diagrams of Figure 1: in the case of synergy (Figure 1a), the two sources $Y_i$ and $Y_k$ taken together contribute to the target $Y_j$ with more information than the sum of their individual contributions $({\mathcal{T}}_{ik\to j}>{\mathcal{T}}_{i\to j}+{\mathcal{T}}_{k\to j})$, and the ITE is positive; in the case of redundancy (Figure 1b), the sum of the information amounts transferred individually from each source to the target is higher than the joint information transfer $({\mathcal{T}}_{i\to j}+{\mathcal{T}}_{k\to j}>{\mathcal{T}}_{ik\to j})$, so that the ITE is negative.

2.2. Partial Information Decomposition

An alternative expansion of the joint TE is that provided by the so-called partial information decomposition (PID) [25]. The PID evidences four distinct quantities measuring the unique information transferred from each individual source to the target, measured by the unique TEs ${\mathcal{U}}_{i\to j}$ and ${\mathcal{U}}_{k\to j}$, and the redundant and synergistic information transferred from the two sources to the target, measured by the redundant TE ${\mathcal{R}}_{ik\to j}$ and the synergistic TE ${\mathcal{S}}_{ik\to j}$. These four measures are related to each other and to the joint and individual TEs by the following equations (see also Figure 1c):

$$\begin{array}{c}\hfill {\mathcal{T}}_{ik\to j}={\mathcal{U}}_{i\to j}+{\mathcal{U}}_{k\to j}+{\mathcal{R}}_{ik\to j}+{\mathcal{S}}_{ik\to j},\end{array}$$

(5a)

$$\begin{array}{c}\hfill {\mathcal{T}}_{i\to j}={\mathcal{U}}_{i\to j}+{\mathcal{R}}_{ik\to j},\end{array}$$

(5b)

$$\begin{array}{c}\hfill {\mathcal{T}}_{k\to j}={\mathcal{U}}_{k\to j}+{\mathcal{R}}_{ik\to j}.\end{array}$$

(5c)

In the PID defined above, the terms ${\mathcal{U}}_{i\to j}$ and ${\mathcal{U}}_{k\to j}$ quantify the parts of the information transferred to the target process $Y_j$, which are unique to the source processes $Y_i$ and $Y_k$, respectively, thus reflecting contributions to the predictability of the target that can be obtained from one of the sources alone, but not from the other source alone. Each of these unique contributions sums up with the redundant transfer ${\mathcal{R}}_{ik\to j}$ to yield the information transfer from one source to the target as is known from the classic Shannon information theory. Then, the term ${\mathcal{S}}_{ik\to j}$ refers to the synergy between the two sources while they transfer information to the target, intended as the information that is uniquely obtained taking the two sources $Y_i$ and $Y_k$ together, but not considering them alone. Compared to the IID defined in (3), the PID (5) has the advantage that it provides distinct non-negative measures of redundancy and synergy, thereby accounting for the possibility that redundancy and synergy may coexist as separate elements of information modification. Interestingly, the IID and PID defined in Equations (3) and (5) are related to each other in a way such that:

$${\mathcal{I}}_{ik\to j}={\mathcal{S}}_{ik\to j}-{\mathcal{R}}_{ik\to j},$$

(6)

thus showing that the interaction TE is actually a measure of the ‘net’ synergy manifested in the transfer of information from the two sources to the target.

An issue with the PID (5) is that its constituent measures cannot be obtained through classic information theory simply subtracting conditional MI terms as done for the IID; an additional ingredient to the theory is needed to get a fourth defining equation to be added to (5) for providing an unambiguous definition of ${\mathcal{U}}_{i\to j}$, ${\mathcal{U}}_{k\to j}$, ${\mathcal{R}}_{ik\to j}$ and ${\mathcal{S}}_{ik\to j}$. While several PID definitions have been proposed arising from different conceptual definitions of redundancy and synergy [26,27,29], here, we make reference to the so-called minimum MI (MMI) PID [36]. According to the MMI PID, redundancy is defined as the minimum of the information provided by each individual source to the target. In terms of information transfer measured by the TE, this leads to the following definition of the redundant TE:

$${\mathcal{R}}_{ik\to j}=min\{{\mathcal{T}}_{i\to j},{\mathcal{T}}_{k\to j}\}.$$

(7)

This choice satisfies the desirable property that the redundant TE is independent of the correlation between the source processes. Moreover, it has been shown that, if the observed processes have a joint Gaussian distribution, all previously-proposed PID formulations reduce to the MMI PID [36].

3. Multiscale Information Transfer Decomposition

3.1. Multiscale Representation of Multivariate Gaussian Processes

In the linear signal processing framework, the M-dimensional vector stochastic process ${\mathit{Y}}_n={[Y_{1,n}\cdots Y_{M,n}]}^T$ is classically described using a vector autoregressive (VAR) model of order p:

$${\mathit{Y}}_n=\sum _{k=1}^p{\mathbf{A}}_k{\mathit{Y}}_{n-k}+{\mathit{U}}_n$$

(8)

where $A_k$ are $M\times M$ matrices of coefficients and ${\mathit{U}}_n={[U_{1,n}\cdots U_{M,n}]}^T$ is a vector of M zero mean Gaussian processes with covariance matrix $\mathbf{\Sigma }$≡$\mathbb{E}\left[{\mathit{U}}_n{\mathit{U}}_n^T\right]$ ($\mathbb{E}$ is the expectation operator). To study the observed process $\mathit{Y}$ at the temporal scale identified by the scale factor $\tau$, we apply the following transformation to each constituent process $Y_m,m=1,\dots ,M$:

$${\overline{Y}}_{m,n}=\sum _{l=0}^q{b}_l{Y}_{m,n\tau -l}.$$

(9)

This rescaling operation corresponds to transforming the original process $\mathit{Y}$ through a two-step procedure that consists of the following filtering and downsampling steps, yielding respectively the processes $\tilde{\mathit{Y}}$ and $\overline{\mathit{Y}}$:

$$\begin{array}{cc}\hfill {\tilde{\mathit{Y}}}_n& =\sum _{l=0}^q{b}_l{\mathit{Y}}_{n-l},\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\overline{\mathit{Y}}}_n& ={\tilde{\mathit{Y}}}_{n\tau },n=1,\dots ,N/\tau \hfill \end{array}$$

(10b)

The change of scale in (9) generalizes the averaging procedure originally proposed in [42], which sets $q=\tau -1$ and $b_l=1/\tau$ and, thus, realizes the step of filtering through the simple procedure of averaging $\tau$ subsequent samples. To improve the elimination of the fast temporal scales, in this study, we follow the idea of [43], in which a more appropriate low pass filter than averaging is employed. Here, we identify the $b_l$ as the coefficients of a linear finite impulse response (FIR) low pass filter of order q; the FIR filter is designed using the classic window method with the Hamming window [51], setting the cutoff frequency at $f_{\tau }=1/2\tau$ in order to avoid aliasing in the subsequent downsampling step. Substituting (8) in (10a), the filtering step leads to the process representation:

$${\tilde{\mathit{Y}}}_n=\sum _{k=1}^p{\mathbf{A}}_k{\tilde{\mathit{Y}}}_{n-k}+\sum _{l=0}^q{\mathbf{B}}_l{\mathit{U}}_{n-l}$$

(11)

where ${\mathbf{B}}_l=b_l{\mathbf{I}}_M$ (${\mathbf{I}}_M$ is the $M\times M$ identity matrix). Hence, the change of scale introduces a moving average (MA) component of order q in the original VAR$(p)$ process, transforming it into a VARMA$(p,q)$ process. As we will show in the next section, the downsampling step (10b) keeps the VARMA representation, altering the model parameters.

3.2. State Space Processes

3.2.1. Formulation of State Space Models

State space models are models that make use of state variables to describe a system by a set of first-order difference equations, rather than by one or more high-order difference equations [52,53]. The general linear state space (SS) model describing an observed vector process $\mathit{Y}$ has the form:

$$\begin{array}{cc}\hfill {\mathit{X}}_{n+1}& =\mathbf{A}{\mathit{X}}_n+{\mathit{W}}_n\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill {\mathit{Y}}_n& =\mathbf{C}{\mathit{X}}_n+{\mathit{V}}_n\hfill \end{array}$$

(12b)

where the state Equation (12a) describes the update of the L-dimensional state (unobserved) process through the $L\times L$ matrix $\mathbf{A}$, and the observation Equation (12b) describes the instantaneous mapping from the state to the observed process through the $M\times L$ matrix $\mathbf{C}$. ${\mathit{W}}_n$ and ${\mathit{V}}_n$ are zero-mean white noise processes with covariances $\mathbf{Q}$≡$\mathbb{E}\left[{\mathit{W}}_n{\mathit{W}}_n^T\right]$ and $\mathbf{R}$≡$\mathbb{E}\left[{\mathit{V}}_n{\mathit{V}}_n^T\right]$ and cross-covariance $\mathbf{S}$≡$\mathbb{E}\left[{\mathit{W}}_n{\mathit{V}}_n^T\right]$. Thus, the parameters of the SS model (12) are ($\mathbf{A},\mathbf{C},\mathbf{Q},\mathbf{R},\mathbf{S}$).

Another possible SS representation is that evidencing the innovations ${\mathit{E}}_n={\mathit{Y}}_n-\mathbb{E}\left[{\mathit{Y}}_n\right|{\mathit{Y}}_n^{-}]$, i.e., the residuals of the linear regression of ${\mathit{Y}}_n$ on its infinite past ${\mathit{Y}}_n^{-}={[{\mathit{Y}}_{n-1}^T{\mathit{Y}}_{n-2}^T\cdots ]}^T$ [53]. This new SS representation, usually referred to as the “innovations form” SS model (ISS), is characterized by the state process ${\mathit{Z}}_n=\mathbb{E}\left[{\mathit{X}}_n\right|{\mathit{Y}}_n^{-}]$ and by the $L\times M$ Kalman gain matrix $\mathbf{K}$:

$$\begin{array}{cc}\hfill {\mathit{Z}}_{n+1}& =\mathbf{A}{\mathit{Z}}_n+\mathbf{K}{\mathit{E}}_n\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill {\mathit{Y}}_n& =\mathbf{C}{\mathit{Z}}_n+{\mathit{E}}_n\hfill \end{array}$$

(13b)

The parameters of the ISS model (13) are ($\mathbf{A},\mathbf{C},\mathbf{K},\mathbf{V}$), where $\mathbf{V}$ is the covariance of the innovations, $\mathbf{V}$≡$\mathbb{E}\left[{\mathit{E}}_n{\mathit{E}}_n^T\right]$. Note that the ISS (13) is a special case of (12) in which ${\mathit{W}}_n=\mathbf{K}{\mathit{E}}_n$ and ${\mathit{V}}_n={\mathit{E}}_n$, so that $\mathbf{Q}=\mathbf{K}\mathbf{V}{\mathbf{K}}^T$, $\mathbf{R}=\mathbf{V}$ and $\mathit{S}=\mathbf{K}\mathbf{V}$.

Given an SS model in the form (12), the corresponding ISS model (13) can be identified by solving a so-called discrete algebraic Riccati equation (DARE) formulated in terms of the state error variance matrix $\mathbf{P}$ [45]:

$$\begin{array}{cc}\hfill \mathbf{P}& =\mathbf{A}\mathbf{P}{\mathbf{A}}^T+\mathbf{Q}-(\mathbf{A}\mathbf{P}{\mathbf{C}}^T+\mathbf{S}){(\mathbf{C}\mathbf{P}{\mathbf{C}}^T+\mathbf{R})}^{-1}(\mathbf{C}\mathbf{P}{\mathbf{A}}^T+{\mathbf{S}}^T)\hfill \end{array}$$

(14)

Under some assumptions [45], the DARE (14) has a unique stabilizing solution, from which the Kalman gain and innovation covariance can be computed as:

$$\begin{array}{cc}\hfill \mathbf{V}& =\mathbf{C}\mathbf{P}{\mathbf{C}}^T+\mathbf{R}\hfill \\ \hfill \mathbf{K}& =(\mathbf{A}\mathbf{P}{\mathbf{C}}^T+\mathbf{S}){\mathbf{V}}^{-1},\hfill \end{array}$$

(15)

thus completing the transformation from the SS form to the ISS form.

3.2.2. State Space Models of Filtered and Downsampled Linear Processes

Exploiting the close relation between VARMA models and SS models, first we show how to convert the VARMA model (11) into an ISS model in the form of (13) that describes the filtered process $\widetilde{{\mathit{Y}}_n}$. To do this, we exploit Aoki’s method [50] defining the state process ${\tilde{\mathit{Z}}}_n={[{\mathit{Y}}_{n-1}^T\cdots {\mathit{Y}}_{n-p}^T{\mathit{U}}_{n-1}^T\cdots {\mathit{U}}_{n-q}^T]}^T$ that, together with $\widetilde{{\mathit{Y}}_n}$, obeys the state Equation (13) with parameters ($\tilde{\mathbf{A}},\tilde{\mathbf{C}},\tilde{\mathbf{K}},\tilde{\mathbf{V}}$), where:

$$\tilde{\mathbf{A}}=\left[\begin{array}{cccccccc}{\mathbf{A}}_1& \cdots & {\mathbf{A}}_{p-1}& {\mathbf{A}}_p& {\mathbf{B}}_1& \cdots & {\mathbf{B}}_{q-1}& {\mathbf{B}}_q\\ {\mathbf{I}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M& {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathbf{0}}_M& \cdots & {\mathbf{I}}_M& {\mathbf{0}}_M& {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M\\ {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M& {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M\\ {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M& {\mathbf{I}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\mathbf{0}}_M& \cdots & {\mathbf{0}}_M& {\mathbf{0}}_M& {\mathbf{0}}_M& \cdots & {\mathbf{I}}_M& {\mathbf{0}}_M\end{array}\right]$$

$$\tilde{\mathbf{C}}=\left[\begin{array}{cccccc}{\mathbf{A}}_1& \cdots & {\mathbf{A}}_p& {\mathbf{B}}_1& \cdots & {\mathbf{B}}_q\end{array}\right]$$

$$\tilde{\mathbf{K}}={\left[\begin{array}{cccc}{\mathbf{I}}_M& {\mathbf{0}}_{M\times M(p-1)}& {\mathbf{B}}_0^{-T}& {\mathbf{0}}_{M\times M(q-1)}\end{array}\right]}^T$$

and $\tilde{\mathbf{V}}={\mathbf{B}}_0$$\mathbf{\Sigma }$${\mathbf{B}}_0^T$, where $\tilde{\mathbf{V}}$ is the covariance of the innovations ${\tilde{\mathit{E}}}_n={\mathbf{B}}_0{\mathit{U}}_n$.

Now, we turn to show how the downsampled process ${\overline{\mathit{Y}}}_n$ can be represented through an ISS model directly from the ISS formulation of the filtered process ${\tilde{\mathit{Y}}}_n$. To this end, we exploit recent theoretical findings providing the state space form of downsampled signals (Theorem III in [45]). Accordingly, the SS representation of the process downsampled at scale $\tau$, ${\overline{\mathit{Y}}}_n={\tilde{\mathit{Y}}}_{n\tau }$ has parameters ($\overline{\mathbf{A}},\overline{\mathbf{C}},\overline{\mathbf{Q}},\overline{\mathbf{R}},\overline{\mathbf{S}}$), where $\overline{\mathbf{A}}={\tilde{\mathbf{A}}}^{\tau }$, $\overline{\mathbf{C}}=\tilde{\mathbf{C}}$, $\overline{\mathbf{Q}}={\mathbf{Q}}_{\tau }$, $\overline{\mathbf{R}}=\tilde{\mathbf{V}}$ and $\overline{\mathbf{S}}={\mathbf{S}}_{\tau }$, with ${\mathbf{Q}}_{\tau }$ and ${\mathbf{S}}_{\tau }$ given by:

$$\begin{array}{cc}\hfill {\mathbf{S}}_{\tau }& ={\tilde{\mathbf{A}}}^{\tau -1}\tilde{\mathbf{K}}\tilde{\mathbf{V}}\hfill \\ \hfill {\mathbf{Q}}_{\tau }& =\tilde{\mathbf{A}}{\mathbf{Q}}_{\tau -1}{\tilde{\mathbf{A}}}^T+\tilde{\mathbf{K}}\tilde{\mathbf{V}}{\tilde{\mathbf{K}}}^T,\tau \ge 2\hfill \\ \hfill {\mathbf{Q}}_1& =\tilde{\mathbf{K}}\tilde{\mathbf{V}}{\tilde{\mathbf{K}}}^T,\tau =1.\hfill \end{array}$$

(16)

Therefore, the downsampled process has an ISS representation with state process ${\overline{\mathit{Z}}}_n={\tilde{\mathit{Z}}}_{n\tau }$, innovation process ${\overline{\mathit{E}}}_n={\tilde{\mathit{E}}}_{n\tau }$ and parameters ($\overline{\mathbf{A}},\overline{\mathbf{C}},\overline{\mathbf{K}},\overline{\mathbf{V}}$), where $\overline{\mathbf{K}}$ and $\overline{\mathbf{V}}$ are obtained solving the DARE (14) and (15) for the SS model with parameters ($\overline{\mathbf{A}},\overline{\mathbf{C}},\overline{\mathbf{Q}},\overline{\mathbf{R}},\overline{\mathbf{S}}$).

To sum up, the relations and parametric representations of the original process $\mathit{Y}$, the filtered process $\tilde{\mathit{Y}}$ and the downsampled process $\overline{\mathit{Y}}$ are depicted in Figure 2a. The step of low pass filtering (FLT) applied to a VAR(p) process yields a VARMA($p,q$) process (where q is the filter order, and the cutoff frequency is $f_{\tau }=1/2\tau$); this process is equivalent to an ISS process [50]. The subsequent downsampling (DWS) yields a different SS process, which in turn can be converted to the ISS form solving the DARE. Thus, both the filtered process ${\tilde{\mathit{Y}}}_n$ and the downsampled process ${\overline{\mathit{Y}}}_n$ can be represented as ISS processes with parameters ($\tilde{\mathbf{A}},\tilde{\mathbf{C}},\tilde{\mathbf{K}},\tilde{\mathbf{V}}$) and ($\overline{\mathbf{A}},\overline{\mathbf{C}},\overline{\mathbf{K}},\overline{\mathbf{V}}$) which can be derived analytically from the knowledge of the parameters of the original process (${\mathbf{A}}_1,\dots ,{\mathbf{A}}_p,\Sigma$) and of the filter ($q,f_{\tau }$). In the next section, we show how to compute analytically any measure appearing in the information decomposition of a jointly Gaussian multivariate stochastic process starting from its associated ISS model parameters, thus opening the way to the analytical computation of these measures for multiscale (filtered and downsampled) processes.

3.3. Multiscale IID and PID

After introducing the general theory of information decomposition and deriving the multiscale representation of the parameters of a linear VAR model, in this section, we provide expressions for the terms of the IID and PID decompositions of the information transfer valid for multivariate jointly Gaussian processes. The derivations are based on the knowledge that the linear parametric representation of Gaussian processes given in (8) captures all of the entropy differences that define the various information measures [37] and that these entropy differences are related to the partial variances of the present of the target given its past and the past of one or more sources, intended as variances of the prediction errors resulting from linear regression [15,17]. Specifically, let us denote as $E_{j|j,n}=Y_{j,n}-\mathbb{E}\left[Y_{j,n}\right|Y_{j,n}^{-}]$, $E_{j|ij,n}=Y_{j,n}-\mathbb{E}\left[Y_{j,n}\right|Y_{i,n}^{-},Y_{j,n}^{-}]$ the prediction error of a linear regression of $Y_{j,n}$ performed respectively on $Y_{j,n}^{-}$ and $(Y_{j,n}^{-},Y_{i,n}^{-})$ and as ${\lambda }_{j|j}=\mathbb{E}\left[E_{j|j,n}^2\right]$, ${\lambda }_{j|ij}=\mathbb{E}\left[E_{j|ij,n}^2\right]$, the corresponding prediction error variances. Then, the TE from $Y_i$ to $Y_j$ can be expressed as:

$$T_{i\to j}=\frac{1}{2}ln\frac{{\lambda }_{j|j}}{{\lambda }_{j|ij}}.$$

(17)

In a similar way, the joint TE from $(Y_i,Y_k)$ to $Y_j$ can be defined as:

$$T_{ik\to j}=\frac{1}{2}ln\frac{{\lambda }_{j|j}}{{\lambda }_{j|ijk}},$$

(18)

where ${\lambda }_{j|ijk}=\mathbb{E}\left[E_{j|ijk,n}^2\right]$ is the variance of the prediction error of a linear regression of $Y_{j,n}$ on $(Y_{j,n}^{-},Y_{i,n}^{-},Y_{k,n}^{-})$, $E_{j|ijk,n}=Y_{j,n}-\mathbb{E}\left[Y_{j,n}\right|Y_{i,n}^{-},Y_{j,n}^{-},Y_{k,n}^{-}]$. Based on these derivations, one can easily complete the IID decomposition of TE by computing $T_{k\to j}$ as in (17) and deriving the interaction TE from (3) and the PID decomposition, as well by deriving the redundant TE from (7), the synergistic TE from (6) and the unique TEs from (5).

Next, we show how to compute any partial variance from the parameters of an ISS model in the form of (13) [44,45]. The partial variance ${\lambda }_{j|a}$, where the subscript a denotes any combination of indexes $\in \{1,\dots ,M\}$, can be derived from the ISS representation of the innovations of a submodel obtained removing the variables not indexed by a from the observation equation. Specifically, we need to consider the submodel with state Equation (13b) and observation equation:

$$Y_n^{(a)}={\mathbf{C}}^{(a)}{Z}_n+E_n^{(a)},$$

(19)

where the superscript ${}^{(a)}$ denotes the selection of the rows with indices a of a vector or a matrix. It is important to note that the submodels (13a) and (19) are not in innovations form, but are rather an SS model with parameters ($\mathbf{A},{\mathbf{C}}^{(a)},\mathbf{K}\mathbf{V}{\mathbf{K}}^T,\mathbf{V}(a,a),\mathbf{K}\mathbf{V}(:,a)$). This SS model can be converted to an ISS model with innovation covariance ${\mathbf{V}}^{(a)}$ solving the DARE (14) and (15), so that the partial variance ${\lambda }_{j|a}$ is derived as the diagonal element of ${\mathbf{V}}^{(a)}$ corresponding to the position of the target $Y_j$. Thus, with this procedure, it is possible to compute the partial variances needed for the computation of the information measures starting from a set of ISS model parameters; since any VAR process can be represented at scale $\tau$ as an ISS process, the procedure allows computing the IID and PID information decompositions for the rescaled multivariate process (see Figure 2).

It is worth remarking that, while the general formulation of IID and PID decompositions introduced in Section 2 holds for arbitrary processes, the multiscale extension detailed in Section 3 is exact only if the processes have a joint Gaussian distribution. In such a case, the linear VAR representation captures exhaustively the joint variability of the processes, and any nonlinear extension has no additional utility (a formal proof of the fact that a stationary Gaussian VAR process must be linear can be found in [37]). If, on the contrary, non-Gaussian processes are under scrutiny, the linear representation provided in Section 3.1 can still be adopted, but may miss important properties in the dynamics and thus provide only a partial description. Moreover, since the close correspondence between conditional entropies and partial variances reported in this subsection does not hold anymore for non-Gaussian processes, all of the obtained measures should be regarded as indexes of (linear) predictability rather than as information measures.

4. Simulation Experiment

To study the multiscale patterns of information transfer in a controlled setting with known dynamical interactions between time series, we consider a simulation scheme similar to some already used for the assessment of theoretical values of information dynamics [15,17]. Specifically, we analyze the following VAR process of order M = 4:

$$\begin{array}{cc}\hfill Y_{1,n}& =2{\rho }_{1}cos2\pi f_1{Y}_{1,n-1}-{\rho }_1^2{Y}_{1,n-2}+U_{1,n},\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill Y_{2,n}& =2{\rho }_{2}cos2\pi f_2{Y}_{2,n-1}-{\rho }_2^2{Y}_{2,n-2}+c{Y}_{1,n-1}+U_{2,n},\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill Y_{3,n}& =2{\rho }_{3}cos2\pi f_3{Y}_{3,n-1}-{\rho }_3^2{Y}_{3,n-2}+c{Y}_{1,n-1}+U_{3,n},\hfill \end{array}$$

(20c)

$$\begin{array}{cc}\hfill Y_{4,n}& =b{Y}_{2,n-1}+(1-b)Y_{3,n-1}+U_{4,n},\hfill \end{array}$$

(20d)

where ${\mathit{U}}_n={[U_{1,n}\cdots U_{4,n}]}^T$ is a vector of zero mean white Gaussian noises with unit variance and uncorrelated with each other ($\Sigma$= $\mathit{I}$). The parameter design in Equation (20) is chosen to allow autonomous oscillations in the processes $Y_i$, $i=1,\dots ,3$, obtained placing complex-conjugate poles with modulus ${\rho }_i$ and frequency $f_i$ in the complex plane representation of the transfer function of the vector process, as well as causal interactions between the processes at a fixed time lag of one sample and with strength modulated by the parameters b and c (see Figure 3). In this study, we set the coefficients related to self-dependencies to values generating well-defined oscillations in all processes (${\rho }_1={\rho }_2={\rho }_3=0.95$) and letting $Y_1$ fluctuate at slower time scales than $Y_2$ and $Y_3$ ($f_1=0.1,f_2=f_3=0.025$). We consider four configurations of the parameters, chosen to reproduce paradigmatic conditions of interaction between the processes:

(a)

isolation of $Y_1$ and $Y_2$ and unidirectional coupling $Y_3\to Y_4$, obtained setting $b=c=0$;

(b)

common driver effects $Y_2\leftarrow Y_1\to Y_3$ and unidirectional coupling $Y_3\to Y_4$, obtained setting $b=0$ and $c=1$;

(c)

isolation of $Y_1$ and unidirectional couplings $Y_2\to Y_4$ and $Y_3\to Y_4$, obtained setting $b=0.5$ and $c=0$;

(d)

common driver effects $Y_2\leftarrow Y_1\to Y_3$ and unidirectional couplings $Y_2\to Y_4$ and $Y_3\to Y_4$, obtained setting $b=0.5$ and $c=1$.

With this simulation setting, we compute all measures appearing in the IID and PID decompositions of the information transfer, considering $Y_4$ as the target process and $Y_2$ and $Y_3$ as the source processes. The theoretical values of these measures, computed as a function of the time scale using the IID and the PID, are reported in Figure 4. In the simple case of unidirectional coupling $Y_3\to Y_4$ ($b=c=0$, Figure 4a), the joint information transferred from $(Y_2,Y_3)$ to $Y_4$ is exclusively due to the source $Y_3$ without contributions from $Y_2$ and without interaction effects between the sources $({\mathcal{T}}_{23\to 4}={\mathcal{T}}_{3\to 4}={\mathcal{U}}_{3\to 4},{\mathcal{T}}_{2\to 4}={\mathcal{U}}_{2\to 4}=0,{\mathcal{I}}_{23\to 4}={\mathcal{S}}_{23\to 4}={\mathcal{R}}_{23\to 4}=0)$.

When the causal interactions towards $Y_4$ are still due exclusively to $Y_3$, but the two sources $Y_2$ and $Y_3$ share information arriving from $Y_1$ ($b=0,c=1$; Figure 4b), the IID evidences that the joint information transfer coincides again with the transfer from $Y_3$ (${\mathcal{T}}_{23\to 4}={\mathcal{T}}_{3\to 4}$), but a non-trivial amount of information transferred from $Y_2$ to $Y_4$ emerges, which is fully redundant (${\mathcal{T}}_{2\to 4}=-{\mathcal{I}}_{23\to 4}$). The PID highlights that the information from $Y_3$ to $Y_4$ is not all unique, but is in part transferred redundantly with $Y_2$, while the unique transfer from $Y_2$ and the synergistic transfer are negligible.

In the case of two isolated sources equally contributing to the target ($b=0.5,c=0$, Figure 4c), the IID evidences the presence of net synergy and of identical amounts of information transferred to $Y_4$ from $Y_2$ or $Y_3$ (${\mathcal{I}}_{23\to 4}>0,{\mathcal{T}}_{2\to 4}={\mathcal{T}}_{3\to 4}$). The PID documents that there are no unique contributions, so that the two amounts of information transfer from each source to the target coincide with the redundant transfer, and the remaining part of the joint transfer is synergistic (${\mathcal{U}}_{2\to 4}={\mathcal{U}}_{3\to 4}=0$, ${\mathcal{T}}_{2\to 4}={\mathcal{T}}_{3\to 4}={\mathcal{R}}_{23\to 4},{\mathcal{S}}_{23\to 4}={\mathcal{T}}_{23\to 4}-{\mathcal{R}}_{23\to 4}$).

Finally, when the two sources share common information and contribute equally to the target ($b=0.5,c=1$; Figure 4d), we find that they send the same amount of information as before, but in this case, no unique information is sent by any of the sources (${\mathcal{T}}_{2\to 4}={\mathcal{T}}_{3\to 4},{\mathcal{U}}_{2\to 4}={\mathcal{U}}_{3\to 4}=0$). Moreover, the nature of the interaction between the sources is not trivial and is scale dependent: at low time scales, where the dynamics are likely dominated by the fast oscillations of $Y_1$, the IID reveals net redundancy, and the PID shows that the redundant transfer prevails over the synergistic (${\mathcal{I}}_{23\to 4}<0,{\mathcal{R}}_{23\to 4}>{\mathcal{S}}_{23\to 4}$); at higher time scales, where fast dynamics are filtered out and the slow dynamics of $Y_2$ and $Y_3$ prevail, the IID reveals net synergy, and the PID shows that the synergistic transfer prevails over the redundant (${\mathcal{I}}_{23\to 4}>0,{\mathcal{S}}_{23\to 4}>{\mathcal{R}}_{23\to 4}$).

5. Application

As a real data application, we analyze intracranial EEG recordings from a patient with drug-resistant epilepsy measured by an implanted array of $8\times 8$ cortical electrodes and two left hippocampal depth electrodes with six contacts each. The data are available in [54], and further details on the dataset are given in [55]. Data were sampled at 400 Hz and correspond to 10-s segments recorded in the pre-ictal period, just before the seizure onset, and 10 s during the ictal stage of the seizure, for a total of eight seizures. Defining and locating the seizure onset zone, i.e., the specific location in the brain where the synchronous activity of neighboring groups of cells becomes so strong so as to be able to spread its own activity to other distant regions, is an important issue in the study of epilepsy in humans. Here, we focus on the information flow from the sub-cortical regions, probed by depth electrodes, to the brain cortex. In [21], it has been suggested that Contacts 11 and 12, in the second depth electrode, are mostly influencing the cortical activity; accordingly, in this work, we consider Channels 11 and 12 as a pair of source variables for all of the cortical electrodes and decompose the information flowing from them using the multiscale IID and PID here proposed, both in the pre-ictal stage and in the ictal stage. An FIR filter with $q=12$ coefficients is used, and the order p of the VAR model is fixed according to the Bayesian information criterion. In the analyzed dataset, the model order assessed in the pre-ictal phase was $p=14.61\pm 1.07$ (mean ± std. dev.across 64 electrodes and eight seizures) and during the ictal phase decreased significantly to $p=11.09\pm 3.95$.

In Figure 5, we depict the terms of the IID applied from the two sources (Channels $\{11,12\}$) to any of the electrodes as a function of the scale $\tau$, averaged over the eight seizures. We observe a relevant enhancement of the joint TE during the seizure, w.r.t. the pre-ictal period. This enhancement is determined by a marked increase of both the individual TEs from Channels 11 and 12 to all of the cortical electrodes; the patterns of the two TEs are similar to each other in both stages. The pattern of interaction information transfer displays prevalent redundant transfer for low values of $\tau$ and prevalent synergistic transfer for high $\tau$, but the values of the interaction TE have relatively low magnitude and are only slightly different in pre-ictal and ictal conditions. It is worth stressing that at scale $\tau$, the algorithm analyzes oscillations, in the time series, slower than $\frac{1}{2\tau f_s}$ s, where $f_s=400$ Hz.

In Figure 6, we depict, on the other hand, the terms of the PID computed for the same data. This decomposition shows that the increased joint TE across the seizure transition seen in Figure 5a is in large part the result of an increase of both the synergistic and the redundant TE, which are markedly higher during the ictal stage compared with the pre-ictal. This explains why the interaction TE of Figure 5d, which is the difference between two quantities that both increase, is nearly constant moving from the pre-ictal to the ictal stage. The quantity that, instead, clearly differentiates between Channels 11 and 12 is the unique information transfer: indeed, only the unique TE from Channel 12 increases in the ictal stage, while the unique TE from Channel 13 remains at low levels.

In order to investigate the variability across trials of the estimates of the various information measures, in Figure 7, we depict the terms of both IID and PID expressed for each ictal episode as average values over all 64 cortical electrodes. The analysis shows that the higher average values observed in Figure 5 and Figure 6 at Scales 1–4 during the ictal state for the joint TE, the two individual TEs, the redundant and synergistic TEs and the unique TE from depth Channel 12 are the result of an increase of the measures for almost all of the observed seizure episodes.

These findings are largely in agreement with the increasing awareness that epilepsy is a network phenomenon that involves aberrant functional connections across vast parts of the brain on virtually all spatial scales [56,57]. Indeed, our results document that the occurrence of seizures is associated with a relevant increase of the information flowing from the subcortical regions (associated with the depth electrode) to the cortex and that the character of this information flow is mostly redundant both in the pre-ictal and in the ictal state. Here, the need for a multiscale approach is testified by the fact that several quantities in the ictal state (e.g., the joint TE, the synergistic ITand the unique ITfrom Channel 12) attain their maximum at scale $\tau >1$.

Moreover, the approaches that we propose for information decomposition appear useful to improve the localization of epileptogenic areas in patients with drug-resistant epilepsy. Indeed, our analysis suggests that Contact 12 is the closest to the seizure onset zone, and it is driving the cortical oscillations during the ictal stage, as it sends unique information to the cortex. On the other hand, to disentangle this effect, it has been necessary to include also Channel 11 in the analysis and to make the PID of the total information from the pair of depth channels to the cortex; indeed, the redundancy between Channels 11 and 12 confounds the informational pattern unless the PID is performed.

6. Conclusions

Understanding how multiple inputs may combine to create the output of a given target is a fundamental challenge in many fields, in particular in neuroscience. Shannon’s information theory is the most suitable frame to cope with this problem and thus to assess the informational character of multiplets of variables describing complex systems; IID indeed measures the balance between redundant and synergetic interaction within the classical multivariate entropy-based approach. Recently Shannon’s information theory has been extended, in the PID, so as to provide specific measures for the information that several variables convey individually (unique information), redundantly (shared information) or only jointly (synergistic information) about the output.

The contribution of the present work is the proposal of an analytical frame where both IID and PID can be exactly evaluated in a multiscale fashion, for multivariate Gaussian processes, on the basis of simple vector autoregressive identification. In doing this, our work opens the way for both the theoretical analysis and the practical implementation of information modification in processes that exhibit multiscale dynamical structures. The effectiveness of the proposed approach has been demonstrated both on simulated examples and on real publicly-available intracranial EEG data. Our results provide a firm ground to the multiscale evaluation of PID, to be applied in all applications where causal influences coexist at multiple temporal scales.

Future developments of this work include the refinement of the SS model structure to accommodate the description of long-range linear correlations [58] or its expansion to the description of nonstationary processes [59] and the formalization of exact cross-scale computation of information decomposition within and between multivariate processes [60]. A major challenge in the field remains the generalization of this type of analysis to non-Gaussian processes, for which exact analytical solutions or computationally-reliable estimation approaches are still lacking. This constitutes a main direction for further research, because real-world processes display very often non-Gaussian distributions, which would make an extension to nonlinear models or model-free approaches beneficial. The questions that are still open in this respect include the evaluation of proper theoretical definitions of synergy or redundancy for nonlinear processes [25,26,27,28,29], the development of reliable entropy estimators for multivariate variables with different dimensions [6,35,61] and the assessment of the extent to which non-linear model-free methods really outperform the linear model-based approach adopted here and in previous investigations [62].