A New Framework for the Time- and Frequency-Domain Assessment of High-Order Interactions in Networks of Random Processes

While the standard network description of complex systems is based on quantifying the link between pairs of system units, higher-order interactions (HOIs) involving three or more units often play a major role in governing the collective network behavior. This work introduces a new approach to quantify pairwise and HOIs for multivariate rhythmic processes interacting across multiple time scales. We define the so-called O-information rate (OIR) as a new metric to assess HOIs for multivariate time series, and present a framework to decompose the OIR into measures quantifying Granger-causal and instantaneous influences, as well as to expand all measures in the frequency domain. The framework exploits the spectral representation of vector autoregressive and state space models to assess the synergistic and redundant interaction among groups of processes, both in specific bands of interest and in the time domain after whole-band integration. Validation of the framework on simulated networks illustrates how the spectral OIR can highlight redundant and synergistic HOIs emerging at specific frequencies, which cannot be detected using time-domain measures. The applications to physiological networks described by heart period, arterial pressure and respiration variability measured in healthy subjects during a protocol of paced breathing, and to brain networks described by electrocorticographic signals acquired in an animal experiment during anesthesia, document the capability of our approach to identify informational circuits relevant to well-defined cardiovascular oscillations and brain rhythms and related to specific physiological mechanisms involving autonomic control and altered consciousness. The proposed framework allows a hierarchically-organized evaluation of time- and frequency-domain interactions in dynamic networks mapped by multivariate time series, and its high flexibility and scalability make it suitable for the investigation of networks beyond pairwise interactions in neuroscience, physiology and many other fields.


I. INTRODUCTION
T HE increasing availability of large-scale and fine-grained datasets is nowadays boosting the development of new methods for the data-driven modelling of complex systems. Among them, the network representation is probably the most used approach for the description of the multivariate time series measured from these systems [1]. Paradigmatic instances of this approach come, among many other fields, from neuroscience and physiology, where the functional connections among different brain regions or among different organ systems are pervasively investigated in the emerging fields of Network Neuroscience [2] and Network Physiology [3]. In this context, data-driven methods for the inference and analysis of complex networks are based on building a network model out of a set of observed time series, in which nodes represent the units composing the observed system (being, e.g., distinct neural populations or physiologic systems) and connecting edges map functional dependencies between units (descriptive, e.g., of brain connectivity or cardiovascular interactions) [4], [5]. Functional dependencies are typically assessed by computing pairwise measures, i.e. measures that describe interactions between two nodes of the analyzed network, on the time series reflecting the dynamic activity of the nodes. The formulation of these measures stems from the availability of several theoretical approaches which formalize the interaction between variables or processes in a network, including the concepts of multivariate spectral analysis [6], [7], Granger causality [8], [9], [10], [11], and directed information transfer and information flow [12], [13], [14], applied to computational systems in the brain and This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in physiology, but also in other cross-disciplinary fields [15], [16], [17], [18], [19], [20], [21], [22], [23].
Nevertheless, in spite of the ubiquitous utilization of pairwise measures to describe interactions in a network, there is mounting evidence that such measures cannot full capture the interplay among the multiple units of a complex system [24]. In fact, complex networks very often exhibit collective behaviors which are integrated at different hierarchical levels, thus displaying interactions that involve more than two network nodes. These socalled high-order interactions (HOIs) occur for instance when brain dynamics require the joint examination of multiple units to be predicted accurately [25], or when cardiovascular interactions are influenced by the effects of the respiratory activity [15].
The recognized need to study networks beyond the framework of pairwise interactions calls for the theoretical definition and practical development of methods to assess HOIs among multiple time series. Various metrics solidly grounded in the general field of information theory have been proposed in recent years for this purpose, all attempting to capture the redundant or synergistic information shared by groups of random variables or processes [25], [26], [27], [28], [29], [30], [31]. In broad terms, synergy arises from statistical interactions that can be found collectively in a network but not in parts of it considered separately, while redundancy refers to group interactions that can be explained by the communication of sub-groups of variables. The most popular measures of synergy and redundancy are those based on interaction information (II, [26]) and on the partial information decomposition (PID, [27]) of random variables, also extended to assess directed interactions in dynamic physiological processes [25], [29]. The II is the first measure proposed to detect synergy and redundancy through their overall balance (the measure is positive when redundancy prevails over synergy, and negative in the opposite case) [26]; the PID provides a different perspective, returning separate and non-negative measures of synergy and redundancy [27] at the cost of not being uniquely defined and difficult to generalize to more than three variables or processes [23], [30]. A recently-proposed measure is the so-called O-information (OI), a metric which extends the II to reveal synergy-and redundancy-dominated interactions in a network of multiple interacting variables [31]. Its symmetric nature, the fact that it scales nicely with the network size, and the possibility to compute it for dynamic processes make the OI a very promising tool for the practical analysis of multivariate dynamics [32].
A main limitation of the information-theoretic measures proposed so far to investigate HOIs in network systems is that they characterize the system dynamics with one single value reflecting the aggregate effect of interactions possibly occurring at different time scales. However, the time series measured at the nodes of complex networks are typically rich of oscillatory content: for instance, cardiovascular and electroencephalographic (EEG) interactions occur through the coupling of rhythms in different frequency bands with different physiological meaning [9], [33]. Remarkably, the amplitude of oscillations and the coupling strength may vary with frequency, and HOIs can have different nature for different rhythms because synergistic and redundant behaviors may alternate in separate frequency bands [7], [17].
Therefore, there is the need to connect the spectral representation of information-theoretic measures with the HOI description of complex networks to overcome spectral pairwise approaches [34], [35]. To this end, the present study introduces a new framework for the time-and frequency-domain analysis of HOIs in multivariate stochastic processes mapping the activity of network systems. Building on our recent efforts to compute multivariate information measures in the frequency domain [7], [17], we generalize and extend them in many directions. First, we define a new measure, the O-information rate (OIR), which generalizes the mutual information rate (MIR) of bivariate processes using the same rationale whereby the OI generalizes the mutual information (MI) between random variables. Then, we provide both a causal decomposition and a spectral expansion of the OIR, thereby connecting it with well-known and widely used measures of coupling and Granger causality formulated in the time and frequency domains [35]. Causal and spectral measures are defined from the vector autoregressive (VAR) formulation of multivariate Gaussian stochastic processes [36], in a way such that the spectral integration of each frequency domain measure yields the corresponding time domain measure. Further, to allow their closed-form computation, all measures composing the time-and frequency-domain OIR are implemented exploiting the state space (SS) representation of VAR processes [37].
In this paper, the proposed framework is first illustrated on theoretical examples of simulated VAR processes featuring HOIs of different type and order. Then, it is tested in two practical applications of of brain and physiological networks where HOIs are expected to play a crucial role in governing collective dynamics: beat-to-beat variability series of heart period, arterial pressure and respiration measured during a protocol of paced breathing [15], and multi-electrode invasive EEG signals acquired in an animal experiment of altered consciousness [18]. The time-and frequency-domain measures of bivariate and higher-order interactions provided by the framework are collected in the OIR Matlab toolbox, described in the supplemental material of this article and freely available for download at www.lucafaes.net/OIR.html.

A. Theoretical Background
This preliminary section reviews the basic and advanced concepts of information theory, applied to random variables and random processes, that pose the basis for the new framework developed in Section II-B.

1) Entropy Measures for Random Variables:
The main measures of information theory are entropy and mutual information (MI), which quantify respectively the information contained in a random variable V 1 , and the information shared by two variables V 1 and V 2 , elaborating their probability distributions as follows: where p(•) and p(•, •) denote joint and marginal probabilities, and E[•] is the statistical expectation operator. The two quantities are linked by the known equation While the MI in (2) quantifies the interaction between two variables, the interaction information (II) is a long-known measure quantifying the interaction among three variables [26], comparing the information shared by one variable with the two other variables when the latter are taken individually or when they are taken together: To perform information-theoretic analysis of HOIs, the II has been recently generalized to quantify interactions among an arbitrarily large number of random variables through the introduction of the so-called O-information (OI) [31]. The OI of a group of N random variables, V N = {V 1 , . . . , V N }, is defined elaborating the entropy of subsets of V N as follows: The OI is a symmetric measure assessing the "organization structure" of a group of random variables; it reduces to the II when evaluated for three variables (i.e., Ω(V 3 ) = I(V 1 ; V 2 ; V 3 )).
The II and the OI defined in (3) and (4) are symmetric measures capturing the balance between high-and low-order statistical constraints in the interactions occurring within V N : Ω(V N ) > 0 reflects a dominance of low-order constraints, also known as redundancy, while Ω(V N ) < 0 indicates that highorder constraints prevail, denoting synergy [26], [31].
2) Entropy Measures for Random Processes: The information measures reviewed above suffer from the limitation that they only allow a static analysis of random variables where the temporal information is disregarded. To perform a dynamic analysis one needs to consider random processes, intended as collections of random variables sorted in temporal order. The generic random process X i is composed by the random variables X i (t n ), where n ∈ N is the temporal index; typically t n = nΔt, where Δt = 1/f s , with f s the sampling frequency. To highlight the dynamic nature of the process, we denote as X i (t n ), X i (t n−k:n−1 ), and X i (t <n ) = lim k→∞ X i (t n−k:n−1 ) the random variables that sample the process at the present time n, over the past k lags, and over the whole past history, respectively. Then, under the assumption of stationarity, the information contained in X i is given by the entropy rate, which quantifies the density of the average information in the process as [38]: Moreover, if two processes X i and X j are considered, the information shared by the processes per unit of time is the mutual information rate (MIR) defined as [39] I X i ;X j = lim k→∞ 1 k I(X i (t n−k:n−1 ); X j (t n−k:n−1 )). (6) Note that, with our notation, H(•) and I(•; •) denote the entropy and MI for random variables, while H (.) and I (.;.) denote the entropy and MI rates for random processes. The entropy rate of a process can be formulated as as the conditional entropy of the present of the process given its past, i.e. [38]. Moreover, starting from the fact that the MIR can be formulated in terms of entropy rates as [38], some elaborations (see e.g. [19], [35]) lead to the important expansion and ) are the transfer entropy (TE) from X j to X i and from X i to X j , and I X i .X j = I(X i (t n ); X j (t n )|X j (t <n ), X j (t <n )) represents the instantaneous information shared between X i and X j ; I(•; •|•) denotes conditional MI for three random variables. The TE is a wellknown measure of directed information transfer between two stochastic processes [40], while the instantaneous transfer is a symmetric measure of information shared at zero lag, quantified after removing the common information with the past states of the processes.

B. Framework to Measure High-Order Interactions in Multivariate Processes
This section presents the formulation of the framework developed to measure dynamic interactions among Q stochastic processes Y Q = {Y 1 , . . . , Y Q }, grouped in M blocks X M = {X 1 , . . . , X M } which can be thought as descriptive of the activity of a network formed by M dynamic systems (the ith block has dimension M i , so that Q = M i=1 M i ). With reference to the applications reported in Section IV, the different dynamic systems analyzed may be M brain regions or M organ systems, where each group process X i , i = 1, . . . , M, may represent the neural activity of a given brain region or organ system, and each scalar process Y j ∈ X i , j = 1, . . . , M i , maps the time course of the jth neural signal recorded inside the ith region (e.g., the EEG at one frontal electrode) or the jth physiological time series belonging to the ith organ system (e.g., systolic or diastolic pressure for the circulatory system). In the following subsections, we define the O-information rate (OIR) as a new measure to assess HOIs among processes, elaborate its causal decomposition, implement its computation in the frame of linear parametric models, and provide its spectral expansion. The framework, whose schematic description is depicted in Fig. 1, allows to study pairwise and higher-order interactions among the analyzed processes both in specific frequency bands related to meaningful rhythmic activities (e.g., brain waves or cardiovascular oscillations) or considering the overall dynamics in the time  (7), initializing to zero the OIR for two processes and then implementing a cycle where the OIR of N processes, Ω X N , is computed adding to the OIR of N − 1 processes, Ω X N −1 , the gradient relevant to the addition of the N th process X N , Δ X N ;X N −1 ; the cycle stops when the OIR of the M processes, Ω X M , is obtained. (b) Iterative computation of the Spectral OIR for the M vector processes. The procedure follows the same steps of the time-domain procedure in (a), applied to the spectral OIR functions ν X N −1 (ω) and ν X N (ω); the core of the procedure is the computation of the OIR gradient δ X N ; which is obtained as a linear combination of N of mutual information rate (MIR) functions (20). Importantly, each time-domain measure is obtained as half the integral of the corresponding spectral function over the whole frequency range. (c) computation of the spectral MIR for a given pair of processes Z 1 and Z 2 . After identifying a vector autoregressive model (VAR) from the Q original processes and converting it into a state space model (SS), a submodel is extracted which contains the parameters relevant only to Z = {Z 1 , Z 2 }; the submodel is analyzed in the frequency domain to derive the spectral measures of Granger and causality and instantaneous interaction that compose the spectral MIR f Z 1 ;Z 2 (ω) according to (23).
domain (e.g., related to brain connectivity or cardiorespiratory coupling).
For our analysis, the processes are assumed to be stationary and ergodic, to allow the time-independent computation of dynamic information measures from individual process realizations [12], [29], and jointly Gaussian distributed, to exploit the formalism linking information-theoretic measures with linear regression models [10], [29] and spectral quantities [7], [17], [35].

1) O-Information Rate:
While the MIR defined in (6) is a dynamic measure of pairwise interdependence between two random processes, HOIs can be assessed generalizing to multiple random processes the OI measure defined in (4) for multiple random variables. Here, following recent works [31], [32], we measure the organization structure of a group of stationary stochastic processes introducing the so-called Oinformation rate (OIR). Specifically, the OIR of the analyzed group of M processes, Ω X M , is defined via the recursion (see also Fig. 1(a)) and where the variation of the OIR obtained with the addition of X N to X N −1 is the quantity with While the OIR can be defined as in (4) using entropy rates in place of entropies, the equivalent formulation (7) highlights the possibility of an iterative computation and evidences the OIR gradient (8) which takes a main role in such computation (see Fig. 1(a)). The OIR is a symmetric measure quantifying redundant and synergistic HOIs among the processes in X N respectively when Ω X N > 0 and Ω X N < 0. In turn, the sign of the OIR gradient detects the informational character of the circuits which link the N th process with the remaining N − 1 processes: the information that Note that when N = 3 processes X 3 = {X 1 , X 2 , X 3 } are considered, substituting (7a) in (7b) yields Ω X 3 = Δ X 3 ;{X 1 ,X 2 } , which expanded with (8) gives a dynamic version of the II measure defined in (3), which we denote as interaction information rate (IIR):

2) Causal Decomposition of the O-Information Rate:
To decompose the OIR increment into causal and instantaneous contributions, we note that Δ X N ;X N −1 is obtained inserting N different MIR values in (8), i.e. the MIRs between the processes . Then, using Z 1 and Z 2 in the MIR expansion [35] and substituting into (8) allows to decompose the OIR gradient as where the three terms quantify the informational character of the directed information transfer from X N to X N −1 , of the directed information transfer from X N −1 to X N , and of the instantaneous information shared between X N −1 and X N , respectively; the informational character of each term is redundant when the term is positive, and synergistic when the term is negative.
3) Linear Parametric Formulation: This subsection reports the parametric implementation of the OIR decomposition, which exploits the knowledge that linear regression models capture all of the entropy differences relevant to the various information measures when the observed processes have a joint Gaussian distribution [10], [29]. As a first step, the analyzed set of stochastic processes Y Q is described as a vector autoregressive (VAR) process of order p: where ] is a Q-dimensional vector random variable collecting the present state of all processes, A(k) is the Q × Q matrix of the model coefficients relating the present with the past of the processes assessed at lag k, . While the VAR model (13) provides a global representation of the overall multivariate process, to describe the linear interactions relevant to the subset of processes Z = {Z 1 , Z 2 } = {X N , X N −1 −i } for which the MIR decomposition is sought we need to define a reduced VAR model involving only those processes. This reduced model is formulated as where An issue with great practical relevance is that the order of the reduced model (14) is typically infinite and thus very difficult to identify from finite-length time series. The approach followed to face this issue in the context of Granger causality analysis is essentially based on truncating to p the order of the reduced model, and estimating its parameters from the relevant subset of the original data. Though simple, this approach exposes to a trade-off between bias and variance of the estimates that prevents reliable model identification in most cases [41]. To solve this issue, methods which essentially extract the parameters of the reduced model from those of the full model have been proposed [42], [43]. Along this line, we overcome the issue related to the formation of the reduced models working in the frame of SS models [37]. This class of models is the most appropriate to use because it is closed under the formation of reduced models: in fact, any reduced process obtained from the VAR process (13) is actually a VAR process with a moving average component, or equivalently a finite-order SS process [43]. Therefore, using SS models allows to identify reduced models from the parameters of the original VAR model estimated with a single regression, thus guaranteeing high computational reliability.
Here, we exploit the SS modeling approach to compute all the MIR terms needed to derive the OIR (7) and to perform the related causal decomposition (11), (12) without the need of re-identifying the parameters of the reduced models from subsets of data ( Fig. 1(c)). First, we describe the original process Y obeying the VAR representation (13) using the SS model is the pQdimensional state process and the SS parameters (A, C, K, V) are given by the matrices  (14) with a reduced SS model with state equation (15a) and observation equation Z(t n ) = C (r,:) S(t n ) + W (t n ). This model has parameters (A, C (r,:) , KVK , V (r,r) , KV (:,r) ), where the superscripts denote selection of the rows and/or columns with indices r in a matrix. To exploit the reduced SS model for the Granger-causal analysis of Z it is necessary to lead its form back to that of (15), which reads [37] The parameters of the reduced model (16) are (Ã,C,K,Ṽ), of dimension pQ × pQ, R × pQ, pQ × R, R × R, and can be derived directly from the parameters A(k) and Σ U of the original full VAR model (13) [37]: while the state and observation matrices are easily determined asÃ = A and and C = C (r,:) , the gainK and the reduced innovation covariancẽ V = E[W n W n ] = Σ W must be obtained by solving a discrete algebraic Riccati equation (DARE) (see [23], [37] for detailed derivations). After identification, the model (16) is analyzed in the frequency domain to compute the spectral components of the MIR, as well as their time-domain counterparts through spectral integration, as reported in the next subsection (see also Fig. 1(b)).

4) Frequency Domain Expansion:
The linear parametric representation of the dynamic interactions among the observed processes can be translated in the frequency domain, in order to provide spectral equivalents of the MIR and OIR measures and of their causal decompositions. Starting from the subset Z = {Z 1 , Z 2 } of the observed multivariate process, described by the SS model (16), taking the Fourier Transform (FT) of the state equation (16a) yields where S(ω) and W (ω) are the Fourier transforms of Z(t n ) and (17) it is easy to derive the PSD of the state process, S(ω), to be substituted in the FT of (16b) to obtain Z(ω) = H(ω)W (ω), which evidences the transfer function matrix The R × R matrix H(ω) contains the transfer functions relating the FTs of the innovation processes in W to the FTs of the processes in Z, and can be used together with the innovation covariance matrix to derive the R × R power spectral density (PSD) matrix of the process Z using spectral factorization: The matrix S Z (ω) can be then factorized in blocks to make explicit the power spectral densities of Z 1 and Z 2 , S Z 1 (ω) and S Z 2 (ω), as diagonal blocks, and the cross-spectral densities between Z 1 and Z 2 , S Z 1 Z 2 (ω) and S Z 2 Z 1 (ω), as off-diagonal blocks. From this factorization, a logarithmic spectral measure of the interdependence between Z 1 and Z 2 is defined by [34] f Z 1 ;Z 2 (ω) = log this measure quantifies the total (symmetric) coupling between Z 1 and Z 2 and is related to the so-called block coherence [44]. Moreover, after factorizing in R i × R i diagonal blocks and R i × R j off-diagonal blocks also the transfer and innovation covariance matrices H(ω) and Σ W , logarithmic spectral measures of the causal effect of Z j on Z i (i, j = 1, 2) can be computed as [34] f where H ii describes the transfer from W i to Z i in the frequency domain and ; these measures quantify the causal (asymmetric) coupling from Z 1 to Z 2 and vice-versa, and are related to the so-called block directed coherence [11].
To complete the representation of the pairwise interactions between Z 1 and Z 2 , a spectral measure f Z i .Z j (ω) can be defined subtracting the sum of the two causal measures (21) from the coupling measure (20) to get so as to satisfy in the frequency domain a decomposition similar to the time-domain decomposition (10): Importantly, the spectral measures in (23) are tightly linked to the similar measures given in the time domain in (10). In fact, it can be shown (see, e.g., [35]) that integration over the whole frequency axis of the spectral coupling measure (20) returns, with proper scaling, the MIR between the two processes, i.e.
and that the same relation holds integrating f Z 1 →Z 2 (ω), This spectral integration property gives to the measures f Z 1 ;Z 2 (ω) and f Z j →Z i (ω) the information-theoretic meaning of density of information shared between the two processes, or transferred from one process to the other, at the angular frequency ω. We note that, while the coupling measure is always non-negative, the two causal measures can take negative values at some frequencies if the process Z is not strictly causal (i.e. if the innovation covariance Σ W is not block-diagonal). On the contrary, the measure f Z 1 .Z 2 (ω) can take negative values even for strictly causal processes [45]. The spectral integration property can be exploited not only to compute the time-domain measures in (10) as the integral of the spectral measures in (23), but also to achieve a causal decomposition of the OIR formulated for spectral functions. Indeed, it is easy to show that the frequency-specific OIR increment defined in analogy to (8) as satisfies the spectral integration property, i.e. Δ X N ;X N −1 = (1/4π) π −π δ X N ;X N −1 (ω) dω, and can also be expanded through a causal decomposition similar to (11) as where the three terms on the r.h.s. of (26) are obtained expanding f X N ;X N −1 (ω) and f X N ;X N −1 −i (ω) in (25) according to (23). Moreover, the spectral OIR increment (25) can be used to compute recursively a frequency-domain version of the OIR, in analogy to (7), as (see Fig. 1 which again satisfies the spectral integration property, i.e. Ω X N = (1/4π) π −π ν X N (ω) dω. Therefore, the spectral versions of the HOI measures defined in this section can be meaningfully interpreted as densities of the synergistic/redundant character of the information shared between multiple stochastic processes. To conclude this section it is worth noting that, in the case of N = 3 processes, the spectral OIR (27) is a frequency-domain analogous of the IIR defined in (9), which can be recovered through whole-band integration. This measure has been recently defined for triplets of random processes [17], and also extended to the spectral computation of separate measures of redundancy and synergy within the PID framework [7]. As shown in the theoretical examples of Section III and practical applications of Section IV, the evaluation of the spectral IIR of three processes, and more generally of the spectral OIR of multiple processes, allows to assess the informational character of specific oscillations within circuits of nodes of the analyzed network.

III. THEORETICAL EXAMPLES
In this section, the framework for the computation of pairwise and higher-order interactions in the time and frequency domains is illustrated making use of theoretical examples of simulated multivariate VAR models for which the various measures are computed directly from the known model parameters. These simulations are exploited to show how our measures can be used: (a) to highlight the emergence of patterns of interaction among groups of processes which cannot be traced from pairwise connections; (b) to dissect pairwise and higher-order interactions into causal components which can be related to the topological structure of the underlying network; (c) to ascribe interactions to specific oscillations confined within specific frequency bands; (d) to evidence the presence of circuits dominated by synergy or redundancy, or even by simultaneous synergistic and redundant behaviors coexisting at different frequencies.
Detailed equations and parameter settings are provided in Section II of the supplemental material, alongside with references to the Matlab codes that implement the two simulations.

A. Simulation 1
The first simulation reproduces the trivariate system proposed in [7], adapted to generate realistic cardiovascular and respiratory dynamics. The activity of this system is mapped by a trivariate VAR process defined as in (13) fed by independent Gaussian innovations, for which the parameters are set as illustrated in Fig. 2(a) and explicitly indicated in (S10) of the Supplemental Material. The vector process is studied keeping the three scalar processes separate (M = Q = 3, X = Y ), and assuming sampling frequency f s = 1 (spectral functions are described completely in the frequency range 0-0.5 Hz). The coefficient matrix A is designed to mimic the dynamics of respiration (X 1 ), arterial pressure (X 2 ) and heart period (X 3 ) variability, generating self-dependencies for the processes X 1 and X 2 through the coefficients a 11,k and a 22,k , and imposing causal effects along the directions X 1 → X 2 , X 1 → X 3 and X 2 → X 3 through the coefficients a 21,k , a 31,k and a 32 . Self-dependencies are set to induce oscillations in the respiratory band (∼ 0.35 Hz) for X 1 and in the low-frequency band (∼ 0.1 Hz) for X 1 and particularly for X 2 , while causal effects are set to realize a high-pass filter from X 1 to X 2 , a low-pass filter from X 1 to X 3 and an all-pass configuration from X 2 to X 3 (spectral transfer functions are shown in Fig. 2(a), right); low-and high-pass filtering are achieved through FIR filters of order 20 with cut-off frequency of 0.2 Hz.
The application of our framework to the VAR parameters describing the simulated process leads to the spectral functions depicted in Fig. 2(b), (c). The PSD profiles (Fig. 2(b), diagonal plots) highlight oscillations at ∼ 0.1 Hz and ∼ 0.35 Hz for the three processes. The causal coupling between pairs of processes ( Fig. 2(b), off-diagonal plots) evidences the presence of information flows originating from the first process (nonzero profiles of f X 1 →X 2 , f X 1 →X 3 and f X 2 →X 3 ) and the absence of information flowing back towards it (f X 3 →X 2 = f X 2 →X 1 = f X 3 →X 1 = 0 at each frequency). Note that, given the unidirectional coupling and the absence of instantaneous interactions, in virtue of (23) the three nonzero causal coupling measures are equivalent to the spectral measures of total coupling f X 1 ;X 2 , f X 1 ;X 3 and f X 2 ;X 3 (red curves in Fig. 2(b)); whole-band integration of such measures by (24) leads to the MIR quantifying the total information shared between pairs of processes, whose values result I X 1 ;X 2 = T X 1 →X 2 = 0.28 nats, I X 1 ;X 3 = T X 1 →X 3 = 0.05 nats and I X 2 ;X 3 = T X 2 →X 3 = 0.24 nats. Then, computation of the MIR between one process and the remaining two leads to obtain the OIR via (8), which for this simulation is Ω X 1 ;X 2 ;X 3 = 0.019 nats, denoting a small redundant interaction among the three processes. Importantly, the spectral expansion (Fig. 2(c)) reveals that this small OIR value is the balance between a synergistic interaction at low frequencies (Ω X 1 ;X 2 ;X 3 = −0.15 nats in the band 0.04-0.12 Hz) and a redundant interaction at higher frequencies (Ω X 1 ;X 2 ;X 3 = +0.33 nats in the band 0.31-0.39 Hz). We also highlight that the causal decomposition of the OIR ν X 1 ;X 2 ;X 3 = δ X 1 ;X 2 ,X 3 reveals the unidirectional nature of the OIR increment (i.e., δ X 1 ;X 2 ,X 3 = δ X 1 →X 2 ,X 3 and δ X 2 ,X 3 →X 1 = δ X 1 .X 2 ,X 3 = 0). The opposite OIR values observed in the two frequency bands can be explained by the simulation design (see Fig. 2(a)): synergy and redundancy arise respectively because the flow of information from X 1 to X 3 is entirely mediated by X 2 at the respiratory frequency (the path X 1 → X 3 is blocked by H 31 at ∼ 0.35 Hz), and because such flow occurs via the independent paths X 1 → X 3 and X 2 → X 3 at lower frequencies (the path X 1 → X 2 is blocked by H 21 at ∼ 0.1 Hz).

B. Simulation 2
The second simulation illustrates the possibility offered by our framework to quantify higher-order spectral interactions among multiple blocks of processes whose dynamics resemble those of neurophyiological signals. The simulation extends previous simulations of VAR processes [11], [17] to the analysis of Q = 10 processes organized in M = 5 blocks, with connectivity structure organized as in Fig. 3(a); equations and parameter setting are given in (S11) of the Supplemental Material. The network is designed to simulate three autonomous vector processes X 1 , X 2 and X 3 which generate, through their own subnetwork interactions, a stochastic oscillation resembling the brain α rhythm (∼ 10 Hz) which is transmitted to the central node X 4 ; such node is a sink for the α waves but also acts as a source of oscillatory activity in the β band (∼ 25 Hz), which is transmitted back to X 1 through the passive block X 5 . The presence of the two simulated rhythms and their transmission through the network is documented by the power spectra S X i and by the pairwise coupling measures f X i ;X j reported respectively in red and gray in Fig. 3(b); integration of the coupling measures leads to detect significant MIR values between each pair of processes except X 2 and X 3 .
The analysis of higher-order interactions was performed computing the spectral OIR for all multiplets of order N = 3, 4, 5 (Fig. 3(c)) as well as the corresponding time-domain OIR values obtained integrating the spectral measures over all frequencies or within the α (8-12 Hz) or β (18-30 Hz) bands (Fig. 3(d)). This analysis allows to evidence patterns of interaction which cannot be inferred from lower-order pairwise links. In particular, the presence of independent sources sending information to a common target originates synergistic modes of interaction characterized by negative profiles of the OIR; this is the case for the multiplets including two or three of the source processes X 1 , X 2 , X 3 and one between X 4 and X 5 (e.g., ν X 1 ,X 2 ,X 4 and ν X 1 ,X 2 ,X 3 ,X 4 , red and violet negative OIRs in Fig. 3(c)). On the contrary, chains of interactions including three or more block processes determine redundant modes of dependence characterized by positive OIR values; this occurs when one or two of the sources X 1 , X 2 , X 3 and both the driven processes X 4 and X 5 are included in the analyzed multiplet (e.g., ν X 1 ,X 4 ,X 5 and ν X 1 ,X 2 ,X 4 ,X 5 , green and cyan positive OIRs in Fig. 3(c)). We note also that the OIR is uniformly null for the triplet with independent processes {X 1 , X 2 , X 3 } (gray line in Fig. 3(c), left panel). The computation of the time-domain OIR puts in evidence the purely synergistic or redundant nature of the interactions occurring within the multiplets of order 3 and 4, as documented in Fig. 3(d) by the clearly negative or positive values of the OIRs. Interestingly, the integration within a specific frequency band (α or β) leads to infer which is the rhythm mostly associated with the interactions, which in this simulation occur dominantly in the α band for the synergistic modes with negative OIR, and in both bands with prevalence of β for the redundant modes with positive OIR.
The analysis of the highest-order multiplet incorporating all processes puts clearly in evidence that synergy and redundancy are related to the simulated α and β rhythms, respectively. Indeed, the spectral OIR ν X 5 displays a negative peak at ∼ 10 Hz and a positive peak at ∼ 25 Hz (Fig. 3(c), right panel), and the integration of this spectral function within the α and β bands evidences clearly negative and positive values (grey bars at the right of Fig. 3(d)). This mode is an example of how the coexistence of synergy and redundancy at different frequencies may mask their time domain detection, as in this case the whole-band integration of the spectral OIR gives small negative values which could be difficult to assess in practice.

IV. APPLICATION TO PHYSIOLOGICAL NETWORKS
This section reports the application of the framework for the analysis of multivariate interactions in the time and frequency domain to two different physiological networks, i.e. cardiovascular and respiratory interactions during paced breathing, and neural interactions from ECoG signals in the anesthetized macaque monkey. Full details about the analyzed datasets and additional results are provided in Section III of the supplemental material.

A. Cardiovascular and Respiratory Interactions During Paced Breathing
The analyzed dataset refers to beat-to-beat variability series of respiration (RESP, process X 1 ), systolic arterial pressure (SAP, process X 2 ) and heart period (HP, process X 3 ), synchronously measured in a group of 18 young healthy subjects monitored in the resting supine position during an experimental protocol consisting of four phases: spontaneous breathing (SB) and controlled breathing at 10, 15, and 20 breaths/minute (CB10, CB15, CB20) [15]. The HP, SAP and RESP time series were extracted respectively from the electrocardiogram, noninvasive arterial blood pressure and nasal respiration flow as the sequences of the duration of the cardiac cycle (R-R interval), of the local maximum of the blood pressure signal within each detected cardiac cycle, and of the value of the respiration signal sampled at the onset of each cardiac cycle. This measurement convention implies that instantaneous influences can be described as causal effects from RESP to SAP and HP and from SAP to HP The analysis was performed on stationary segments of the time series including 256 heartbeats, selected by visual inspection for each subject and experimental condition [15]. The pre-processing consisted on detrending and mean removal for each time series. The VAR model fitting the three series was identified through the ordinary least squares method, selecting the order p in the range 3-14 by means of the Akaike Information Criterion [36]. The analysis was focused on decomposing the OIR of the three processes in OIR increments obtained when the HP process is added to the bivariate process {RESP,SAP}. Specifically, starting from the estimated VAR parameters, we computed δ X 1 ,X 2 →X 3 (f ), δ X 3 →X 1 ,X 2 (f ) and δ X 1 ,X 2 .X 3 (f ) from the terms of the spectral decomposition (10), then deriving ν X 1 ,X 2 ,X 3 (f ) = δ X 1 ,X 2 ;X 3 (f ) via (13,14). From these spectral measures, time-domain measures were obtained through integration over the whole frequency axis or within the low frequency range (LF, 0.04-0.12 Hz) and the high frequency range (HF, ±0.04 Hz around the respiratory frequency f RESP ). Given the possibility to ascribe instantaneous effects to specific causal directions (see above), the analysis is performed summing the information shared instantaneously between {RESP,SAP} and HP to the information transferred from { RESP,SAP} to HP, i.e. computing the spectral and time domain measures δ X 1 ,X 2 .
The results of OIR computation and decomposition are reported in Fig. 4, showing the grand average of the frequencydomain measures as well as the whole-band, LF and HF time-domain average measures. Spectral analysis was performed assuming the series as uniformly sampled with sampling frequency equal to the inverse of the mean HP. The spectral OIR and most of the terms of its decomposition exhibit prominent peaks, which are well-defined at the frequency of the paced breathing during the CB conditions and are less narrow-banded during SB (Fig. 4(a)). This behavior reflects the fact that paced breathing regularizes the RESP signal around the imposed rhythm and enforces synchronous oscillations at the same frequency in the HP and SAP time series, determining increased spectral content and spectral coupling in the HF band [15]. The positive values of the time-domain OIR ( Fig. 4(b), left) document that this synchronized interaction is dominantly redundant, confirming previous findings [29].
Looking at the spectral profiles of Fig. 4(a), the peak values of the OIR show a tendency to increase while moving from SB to CB10, and to decrease progressively during CB15 and CB20; these trends confirm from the perspective of HOIs results obtained on the same data using information-theoretic measures of cardiorespiratory coupling [20]. The dominance of redundancy in the HF band of the spectrum (Fig. 4(b), right) suggests that the main underlying physiological mechanism is the mechanical influence of RESP on SAP variability, transmitted to HP through the baroreflex feedback [47]; the OIR component directed from HP to {SAP,RESP}, which tends to be less redundant at increasing the frequency of paced breathing, is of more difficult interpretation and is likely dominated by the mechanical feedforward effects from HP to SAP [21]. The dominance of redundant mechanisms around the respiratory frequency impacts substantially the whole-band time-domain OIR, which show comparable values across the analyzed conditions ( Fig. 4(b), left). On the other hand, the measures integrated within the LF band vary significantly moving from spontaneous to paced breathing ( Fig. 4(b), middle): the information transfer from {SAP,RESP} to HP becomes mostly synergistic during CB10, and during CB15 and CB20 returns progressively to the redundant values observed at SB; the information transfer along the direction HP→ {SAP,RESP} is prevalently synergistic at rest and shifts to redundant values during CB. The shift to synergy observed at CB10 for Δ X 1 ,X 2 →X 3 suggests that, when the respiratory activity slows down and tends to overlap with the Mayer waves typically observed in SAP and HP [48], the baroreflex (SAP→ HP) and respiratory sinus arrhythmia (RESP → HP) mechanisms operate independently in determining the variability of heart rate.

B. Neural Interactions From ECoG Signals in the Anesthetized Macaque Monkey
The second practical application refers to monkey electrocorticographic (ECoG) signals downloaded from the public server neurotycho.org. The analyzed dataset was recorded with a sampling frequency of 1000 Hz in one macaque monkey using 128 electrodes, placed in pairs with an inter-electrode distance of 5 mm to cover the frontal, parietal, temporal and occipital lobes of the left hemisphere [18]. Specifically, we considered two five-minutes recording sessions during which the blindfolded monkey was seated in a primate chair with tied hands, first in a resting state (REST) and then after injection of a sedative inducing anesthesia (ANES). From the 128 electrodes, a subset of 20 was selected as depicted in Fig. 5(a) to cover, considering ten bipolar ECoG signals obtained taking the differential activity between close electrodes, the following five brain regions of the default mode network, i.e. the pre-frontal cortex ( , and high visual cortex (X 5 = [Y 9 , Y 10 ]). The ten bipolar signals were band-pass filtered between 0.5 and 200 Hz, downsampled to f s = 250 Hz, epoched to extract ∼ 160 trials lasting 2 sec for each condition, and finally normalized to zero mean and unit variance within each trial. Then, a VAR model was fitted on the Q = 10 signals of each trial using least squares identification and setting the model order according to the Bayesian Information Criterion (BIC) [36]. From the VAR parameters, the analysis of high-order interactions was performed for the M = 5 blocks computing the spectral OIR for all multiplets of order N = 3, 4, 5. Time-domain OIR values (Ω) were then obtained integrating the spectral measures ν(f ) within the δ (0.2-3 Hz), θ (4-7 Hz), α (8-12 Hz), β (12-30 Hz) and γ (31-70 Hz) frequency bands, as well as cumulatively between 0 and 70 Hz.
The results of OIR computation are reported in Fig. 5(b), showing the grand average of the spectral OIR for five multiplets selected as the most representative of the analyzed interactions, together with the time-domain OIR obtained through wholeband and band-specific integration. The positive values of the OIR functions and of the integrated measures, observed for all multiplets in both conditions and increasing with the order of the multiplet, indicate that the analyzed system is dominated by redundancy. Moreover, the redundancy level is modulated by the experimental condition to an extent that depends on the analyzed multiplet and spectral band. Indeed, considering the multiplets of order 3 and 4 which involve the prefrontal cortex X 1 (1st and 3rd row of panels in Fig. 5(b)), a significant increase of the OIR is observed while moving from REST to ANES; such increase is driven by the rise of a peak in the OIR at ∼ 2 Hz (δ band) together with an increased contribution within the γ band. On the other hand, the multiplets formed by signals from the parietal, temporal and visual cortices (2nd and 4th row of panels in Fig. 5(b)) display a drop of redundancy in the α and β bands during ANES compared to REST. These two opposite behaviors are summarized by the OIR encompassing all five regions (5th row of panels in Fig. 5(b)), which during ANES displays significantly higher levels of redundancy in the δ and γ bands (and in the whole band), and significantly lower redundancy in the θ, α, and β bands.
Our results indicate that the activity relevant to the α and β rhythms observed during the relaxed awake state disappears during anesthesia, leaving place to dominant interactions within the δ and γ bands. The redundancy observed at REST for the α waves is significant for the multiplets involving signals from the visual cortex, in agreement with the knowledge that these waves can be predominantly recorded from the occipital lobes during wakeful relaxation with closed eyes [49]. On the other hand, the higher redundancy reported in the δ band can be related to the slow wave oscillations (0.1-4 Hz) typically observed under anesthesia [50]. Moreover, the fact that higher δ redundancy is observed only for multiplets including frontal cortex signals supports the knowledge that the slow oscillations are a manifestation of a coupling between the anterior and posterior axes of the brain [51]. Anesthesia evokes also an increase of redundancy related to γ oscillations, which are associated with different cognitive functions [52].
Overall, these results agree with those in [18] and support the integration theory according to which the conscious state is generated by highly integrated neural interactions that disappear in the unconscious state [53]. A recent study comparing resting wakefulness with propofol-induced anaesthesia in human fMRI data has shown how the anterior-posterior disconnection occurring during anesthesia is associated with a decrease of Integrated Information within the default mode network in the left hemisphere [54]. Importantly, the concepts of Integration Information and that of redundancy are interrelated, as explained in [55] where it is highlighted that a drop of Integrated Information corresponds to an increase of redundancy. Thus, our results support the theory of an anterior-posterior disconnection during anesthesia, which in our case can be ascribed to the significant increase of the OIR documented when the frontal cortex is considered in the analyzed multiplet.

V. CONCLUSION
This work opens the way to the combined informationtheoretic and spectral evaluation of hierarchically-organized interactions in dynamic networks mapped by multivariate stochastic processes. The proposed framework is highly flexible and scalable as it provides principled measures of both pairwise and higher-order interactions among scalar or vector processes, defined in both time and frequency domains in a way such that the two representations are connected in a straightforward way. Moreover, it allows to decompose symmetric measures into components reflecting Granger-causal and instantaneous influences, and to estimate them with high computational reliability within the framework of vector autoregressive and state space (SS) models.
The application of the new framework to biomedical time series illustrates its capability to capture the balance between redundancies and synergies among arbitrarily large groups of nodes of brain and physiological networks. Moreover, it highlights the importance of studying these features within specific frequency bands of biological interest to elicit interactions which may be otherwise hidden if investigated only in the time domain. The generality of the information-theoretic grounds and of the parametric implementation of the proposed approach makes it suitable for the assessment of pairwise and higher-order interactions even beyond the domain of biomedical time series, to analyze virtually any type of dynamic network (e.g., electronic, climatologic, social, or financial) with node activity described by rhythmic processes.