Measuring High-Order Interactions in Rhythmic Processes through Multivariate Spectral Information Decomposition

Many complex systems in physics, biology and engineering are modeled as dynamical networks and described using multivariate time series analysis. Recent developments have shown that the emergent dynamics of a network system are significantly affected by interactions involving multiple network nodes which cannot be described using pairwise links. While these higher-order interactions can be probed using information-theoretic measures, a rigorous framework to describe them in the frequency domain is still lacking. This work presents an approach for the spectral decomposition of multivariate information measures, capable of identifying higher-order synergistic and redundant interactions between oscillatory processes. We show theoretically that synergy and redundancy can coexist at different frequencies among the output signals of a network system and can be detected only using the proposed spectral method. To demonstrate the broad applicability of the framework, we provide parametric and non-parametric data-efficient estimators for the spectral information measures, and employ them to describe multivariate interactions in three complex systems producing rich oscillatory dynamics, namely the human brain, a ring of electronic oscillators, and the global climate system. In these systems, we show that the use of our framework for the spectral decomposition of information measures reveals multivariate and higher-order interactions not detectable in the time domain. Our results are exemplary of how the frequency-specific analysis of multivariate dynamics can aid the implementation of assessment and control strategies in realworld network systems.


I. INTRODUCTION
The complexity of many physical, biological and technological systems originates from the richness of the structural and The associate editor coordinating the review of this manuscript and approving it for publication was Sung Chan Jun .functional interactions among their constituent units.Science has nowadays abandoned the reductionist notion that the collective behavior of a network system, i.e. a system composed by several possibly interacting units, can be understood and predicted by considering each individual unit in isolation.Moreover, even the classical representation of a network system based on quantifying the interactions between unit pairs through pairwise measures of coupling or causality [1], [2] is being debated, as it is deemed as insufficient to characterize exhaustively the emergent collective dynamics of the system.In fact, as the exploration of real-world systems intensifies, multidisciplinary scientists realize the need to go beyond the framework of pairwise interactions [3].There is indeed mounting evidence that the overall interplay among the several units of a network system cannot be exhaustively described by combinations of pairwise couplings, and that higher-order interactions −i.e., interactions involving more than two units− are present and often play a crucial role for understanding the overall system behavior [4], [5].In this context, novel approaches are under development to accommodate higher-order many body interactions into generalized network representations [6].Nevertheless, the implementation of these generalized network structures remains difficult in systems where interactions are not already identified but need to be inferred from data, as is the case for biological networks and functional networks mapping the dynamics of physical systems.An approach increasingly employed to achieve a proper description of dynamic networks is the use of information theory.
Indeed, information-theoretic measures are being exploited extensively to characterize multivariate dynamical systems across very diverse areas such as neuroscience, finance, physiology and climatology [7]- [9].The framework of information dynamics provides a unifying set of measures which allow quantifying the amounts of information produced and stored in a complex system, as well as those transferred from a ''source'' to a ''target'' and modified as a consequence of the interaction between multiple sources sending information to the target [10], [11].In particular, information modification is related to the concepts of redundancy and synergy between two source systems sharing information about a target system.These concepts address the existence of shared information about the target that can be recovered when the sources are used separately (redundancy), or when they are used jointly (synergy) [12]- [14].Synergy and redundancy have been successfully implemented in measures of interaction information [11], implicitly providing a means to quantify higher-order interactions from the dynamics of complex systems.Interaction information measures were employed to study, among others, cardiovascular and cardiorespiratory interactions [15], electroencephalographic (EEG) brain connectivity [13], and global financial markets [7].While these measures are generally defined in a model-free framework on the basis of the probability densities of the various available dynamic variables, it has been shown that for Gaussian processes they are fully dependent on the parameters of a linear Vector Autoregressive (VAR) model [11], [16].
Like the other measures of information dynamics, interaction information evaluates the overall dynamics of multivariate processes and, as such, operates in the time domain.On the other hand, the signals measured at the nodes of real-world networks are often rich of oscillatory content, and therefore naturally lend themselves to spectral representation.The frequency domain analysis of pairwise coupling is usually performed through spectral measures such as coherence, partial coherence, directed coherence, and partial directed coherence [17], [18].However, these measures do not provide frequency-specific information about higher-order interactions, like interaction information, redundancy and synergy.To fill this gap, we recently introduced a spectral method, relying on VAR models studied in the frequency domain, to assess the information shared among triplets of scalar processes [19].The method decomposes across frequencies the information shared between the processes, providing spectral quantities that −when integrated over the whole frequency axis− return well-known time domain measures of linear dependence [20], [21].This approach opens the way to the evaluation of higher-order interactions in the frequency domain, but is limited in the fact that it is designed only for scalar source and target processes (therefore, for the exclusive analysis of three time series), it performs only VAR parametric estimation of the proposed measures, and it is validated experimentally only in the specific context of cardiovascular variability analysis.

A. SCOPE AND STRUCTURE OF THIS STUDY
The present work aims at making the new framework for the spectral quantification of higher-order interactions [19] widely applicable in many contexts of science and engineering.To this end, we first extend the theoretical formulation to fully multivariate processes, so as to allow the analysis of networks composed by more than three nodes.Moreover, we consider a model-free estimator of the spectral and information measures and we compare it with the parametric estimator based on VAR models in simulations of multivariate time series.Importantly, we also show the broad applicability of the framework in different real-world scenarios: brain networks probed through electroencephalographic (EEG) signals recorded during motor execution [22]; a ring of chaotic electronic oscillators producing non-trivial functional networks [23]; and climate system dynamics assessed from global temperatures and carbon emissions [24].Though apparently rather diverse from each other, these network systems are prime examples of dynamical systems producing broadband signals with evident rhythmic behavior.
Moreover, emergent properties can arise in all these systems from non-trivial interactions within the network, producing global behaviors that are rarely explainable focusing on the activity of individual nodes or the coupling between pairs of nodes.Our focus in these applications is assessing the presence of higher-order interactions quantified by redundant or synergistic behaviors among the system units, and to show that these behaviors can be often elicited only looking at interactions in the frequency domain and even be simultaneously present when observed at different frequencies.
The remainder of this paper is organized as follows.In Sec.II, the theoretical underpinnings of the framework are illustrated and explained.In Sec.III, the behavior of the proposed measures is investigated through theoretical examples, and the performances of parametric and non-parametric estimators for these measures are evaluated.In Sec.IV, the practical applicability of the framework is explored by considering the dynamics of three widely different complex systems.Finally, in Sec.V, the main features of the proposed interaction information decomposition are discussed and compared with other known and widely employed measures and analyses.
The interaction measures, theoretical example, simulation study, and physical scenario analysis are collected in the fdIID MATLAB toolbox available at https://github.com/YuriAntonacci/fdIID-toolbox.Codes can be downloaded also from http://www.lucafaes.net/fdID.html.

II. METHODS
In this section, we first introduce the measures used in the literature to assess dependencies between blocks of stochastic processes in the frequency domain, highlighting the link between such measures and information-theoretic measures.These measures are subsequently exploited to define frequency domain and information-theoretic measures of redundancy/synergy.Then, we describe the procedures for estimating all measures from multivariate time series data, as well as to assess their statistical significance.

A. INFORMATION AND SPECTRAL MEASURES OF LINEAR ASSOCIATION BETWEEN MULTIVARIATE GAUSSIAN PROCESSES
This section aims to set the basic measures and relations which establish the link between information-theoretic and spectral measures of coupling in multivariate (vector) processes.Specifically, we focus on the Mutual Information Rate (MIR) and on its relation with frequency domain measures.Linear dependencies between multivariate processes can be expressed with common statistical tools such as correlation or higher moments measures [25].The approach described in the following relies of second-order statistical measures and their spectral equivalents (e.g.coherence, block coherence).
Let us consider two zero-mean stationary multivariate (vector) stochastic processes Let us also define the present of the processes as the vector containing the random variables which sample X and Y at time An information-theoretic measure of the total degree of association between the two processes is MIR, defined as: [20] where the quantity I (X n ; Y n ) represents the mutual information between the random variables which collect the entire history of X and Y up to time n.
In the frequency domain, the multivariate processes can be described on the basis of the power spectral density (PSD) of each constituent scalar process, defined as the Fourier Transform (FT) of the auto-correlation function of the process (e.g., with k representing the time lag), and of the cross-spectral density between two scalar processes, defined as the FT of their cross-correlation function (e.g., , f s sampling frequency).As customary, the spectral densities are collected in the individual square PSD matrices where P YX (ω) = P * XY (ω) ( * stands for Hermitian transpose).The individual PSD matrices P X (ω) and P Y (ω) and the joint PSD matrix P [XY ] (ω) can be exploited to obtain the spectral measure of total interdependence between X and Y defined as [26]: which quantifies the total interdependence between the two processes in the time domain.Importantly, the measure (3) has an information-theoretic interpretation, since it appears in the spectral representation of the MIR defined in (1).Indeed, a long-known result is that, for Gaussian processes, the MIR can be expanded in the frequency domain as: [20] which also shows that F X ;Y = 2I X ;Y .Given this interpretation, and using the natural logarithm, the quantity F X ;Y is measured in natural units (nats), and the spectral quantity f X ;Y (ω) is measured in nats/rad.In closing this section, we note that the methodological approach used in this work follows the extension of the classical coherence function to the multivariate case through the concept of block coherence [27].In fact, the spectral measure (3) and the block coherence C (b) X ;Y can be related via: ω) ; the equation reduces to the spectral coherence [28] if X and Y are scalar processes.This relationship, together with (5) and considering also the link between cross-correlation and coherence, contributes to relate timedomain, frequency-domain and information-theoretic measures of linear association between stochastic processes [29].

B. SPECTRAL INFORMATION DECOMPOSITION FOR MULTIVARIATE PROCESSES
Now, let us consider a multivariate process where Y is taken as the ''target'' and X 1 and X 2 are assumed as ''sources''.We stress that both the target and the source processes can be multivariate, so we study the interactions between ''blocks'' of time-series considered as realizations of these processes.Let M , M Y , M 1 , and M 2 be the sizes of Z , Y , X 1 and The interaction in the frequency domain can be studied by considering the (M ×M )-dimensional spectral density matrix In (6), P X 1 (ω) is an (M 1 × M 1 )-dimensional matrix containing the PSD of the scalar process X 1i as the i th diagonal element and the cross-PSD between X 1i and X 1j as the (i, j) off-diagonal element, and P X 1 X 2 (ω) is an (M 1 × M 2 )dimensional matrix containing the cross-PSD between X 1i and X 2j as the (i, j) element; the same follows intuitively for the matrices P X 2 (ω), P Y (ω), P X 1 Y (ω) and P X 2 Y (ω), having size M 2 ×M 2 , M Y ×M Y , M 1 ×M Y and M 2 ×M Y .Considering the overall spectral matrix P Z (ω) and its constituent blocks, time-and frequency-domain measures of total dependence between X 1 and Y , between X 2 and Y , and between X = [X 1 X 2 ] and Y , can be defined as in (3) and (4): where the integrands correspond respectively to f X 1 ;Y (ω), f X 2 ;Y (ω) and f X 1 X 2 ;Y (ω).Then, based upon the definition of interaction information for random variables [30], we define the following information-theoretic and spectral measures of source interaction: Since interaction information can assume both positive and negative values, the definitions given in (8) and (9) put in evidence redundant and synergistic interactions in the analyzed system.Specifically, when the two sources X 1 and X 2 , taken individually, are more strongly coupled to the target Y than in the case in which they are taken jointly (F X 1 ;Y + F X 2 ;Y > F X 1 X 2 ;Y ), the time-domain source interaction measure is positive (I X 1 ;X 2 ;Y > 0), indicating redundancy.Conversely, when the two sources X 1 and X 2 taken jointly exhibit a stronger coupling with the target Y than in the case in which they are taken individually (F X 1 ;Y + F X 2 ;Y < F X 1 X 2 ;Y ), the time-domain source interaction measure is negative (I X 1 ;X 2 ;Y < 0), indicating synergy.The same relations hold for the frequency-domain interaction measure, and are thus valid at each specific frequency.In addition, the time-and frequency-domain interaction measures satisfy the property of spectral integration, i.e.I X 1 ;X 2 ;Y = 1 2π π −π i X 1 ;X 2 ;Y (ω) dω.Importantly, we also note that the time-and frequency-domain measures of interaction information ( 8) and ( 9) are symmetric under the exchange of target and source processes.

C. POWER SPECTRAL DENSITY (PSD) ESTIMATION
In the linear signal processing framework, the vector process Z = {Z 1 , . . ., Z M } T can be represented as the output of a multivariate linear shift-invariant filter [31]: where T is a vector of M zero-mean input processes denoted as innovations, and H k is the M × M filter impulse response.A very common representation of (10) is the vector autoregressive (VAR) representation: where p is the model order, defining the maximum lag used to quantify the interactions, and A k are the M ×M coefficient matrices wherein each element a ij,k describes the dependence of Z i,n on Z j,n−k (i, j = 1, . . .M ; k = 1, . . ., p).The innovation process U n contains white and uncorrelated noises, so that it is fully described by its zero-lag covariance matrix A spectral representation can be obtained by taking the discrete-time FT of the representations (10) and (11), which yields Z (ω) = H(ω)U (ω) and Z (ω) = A(ω)Z (ω) + U (ω), where Z (ω) and U (ω) are the FTs of Z n and U n and the M × M transfer matrix and coefficient matrix are defined in the frequency domain as: Comparing the two spectral representations, it is easy to show that the coefficient and transfer matrices are linked VOLUME 9, 2021 by: H(ω) = [I − A(ω)] −1 .Then, by using the spectral factorization theorem [32], it is possible to obtain the spectral density matrix P Z (ω) = H(ω) U H * (ω) (where * stands for complex conjugate).This matrix can be partitioned as in (6) to obtain the frequency domain measures of coupling and interaction, and the corresponding time domain measures through integration over the whole frequency axis, as seen in the previous subsection.
The identification procedure of the VAR model ( 11) is typically performed through the multivariate version of the ordinary least-squares (OLS) method [33].In brief, defining the past history truncated at p time-steps as the pM-dimensional vector Z The method estimates the coefficient matrices through the well-known OLS formula, A = z(z p ) T [z p (z p ) T ] −1 ; then, the innovation process is estimated as the residual time-series U = z − Az p , whose covariance matrix U is an estimate of the innovation covariance.Finally, the spectral density matrix is estimated as P Z (ω) = H(ω) U H * (ω), where the transfer matrix is An alternative to the parametric estimation of the PSD is the weighted covariance (WC) method, which represents a non-parametric approach to deriving the PSD through the FT of the sample autocorrelation and cross-correlation functions of the data [34].Specifically, following the notation introduced for the parametric approach, correlation and power spectral density can be defined for the multivariate process Z as R Z (k) = E[Z n Z T n−k ] and R Z (ω) = F{R Z (k)}, respectively.These are matrices of size M × M that contain the correlation between two scalar processes, R Z i Z j (k), and its FT, as i − j elements (i, j = 1, . . ., M ).The WC estimator of the PSD matrix computes the cross-PSD between the two generic processes Z i and Z j as note that the autocorrelation of Z i and the corresponding PSD are obtained when i = j.In (13), τ ≤ N − 1 is the maximum lag for which the correlation is estimated (with N being the number of data samples available), and w is a lag window of width 2τ (w(k Window selection is usually performed by looking at the spectral leakage introduced by the profile of the window [35].In this work, we used a biased estimator for both cross-correlation and autocorrelation functions which guarantees semi-definite sequences and thus does not lead to negative spectral estimates.A biased estimator of the cross-correlation function can be defined as where the latter holds for k = 0, . . ., N − 1; if k = −(N − 1), . . ., −1, the auto-covariance matrix is defined as

D. TESTING THE SIGNIFICANCE OF INTERACTION INFORMATION IN THE TIME AND FREQUENCY DOMAINS
Since interaction information is a measure of coupling between the sources (X 1 , X 2 ) and the target (Y ), it is of interest to establish its zero-level confidence bounds in order to assess the existence of statistically significant synergistic or redundant contributions in the analyzed network.In this work, the significance of interaction information in the time and frequency domains was tested generating sets of surrogate time-series for two multivariate sources X 1 and X 2 , while leaving the multivariate target process Y untouched.Specifically, surrogates of the two sources were obtained through the phase randomization procedure, preserving not only each individual spectrum but also the magnitude of the cross-spectrum between each pair of series belonging to X 1 and X 2 [36], [37].This is accomplished by adding the same random number to the Fourier phases of the two groups of time-series, so as to maintain autocorrelation and linear coupling between the two groups of sources, and to destroy any coupling between each source and the target process.
The procedure was repeated 100 times, estimating for each repetition the interaction information in the time and frequency domains.The value of interaction information computed on the original data was then compared to the surrogate data, taking as thresholds the 2.5 th and 97.5 th percentiles of the surrogate distributions.The original value was deemed as indicating statistically significant redundant information when it was positive and above the superior threshold.On the other hand, statistically significant synergistic information was inferred when the value of interaction information computed on the experimental data was negative and below the inferior threshold.

III. SIMULATION STUDY
This section reports the theoretical design and the practical implementation of a simulation of multiple interacting stochastic processes, which is used as a benchmark to illustrate the properties of the spectral information measures defined in this work.Using this simulation, we show theoretically that the proposed framework can elicit patterns of redundancy and synergy which coexist at different frequencies in the spectral domain but may be hidden if only the integrated measures in the time domain are computed.We also address the estimation of these measures, showing that both the parametric approach based on VAR models and the non-parametric method based on cross-covariance estimation can provide accurate estimates, with performance depending on the setting of the method-specific free parameters.

A. THEORETICAL EXAMPLE
To study the behavior of the presented measures, we design a four-variate VAR process configured to reproduce coexisting redundant and synergistic interactions between source processes sharing information with a target [8], [38].We consider a theoretical example of four Gaussian systems whose associated processes are described by the VAR with model order equal to 3 (VAR(3)) defined by the following equations: where T is a vector of zero-mean white Gaussian noises with unit variance and uncorrelated with each other ( U = I).The parameter design in ( 15) is chosen to allow autonomous oscillations in the processes Z i , i = 2, 3, 4, obtained placing complex-conjugate poles with modulus ρ 2 = ρ 3 = 0.85 and ρ 4 = 0.95 and normalized frequencies f 2 /f s = 0.05, f 3 /f s = 0.05 and f 4 /f s = 0.35 in the complex plane; assuming a sampling frequency f s = 100 Hz, the poles determine oscillations at 5 Hz and 35 Hz.Moreover, interactions between different processes are set to obtain a common driver effect Z 2 ← Z 4 → Z 3 and unidirectional couplings Z 3 → Z 1 and Z 2 → Z 1 , with weights indicated in Fig. 1(a).Given (15), we take X 1 = Z 2 and X 2 = [Z 3 , Z 4 ] as source processes, and Y = Z 1 as target process.The VAR(3) model described above can be studied in the frequency domain by deriving the 4 × 4 spectral density matrix P Z (ω) directly from the AR coefficients as described in the previous section.This leads to computing the exact values of all the timeand frequency-domain information measures for the theoretical process.Varying the parameter c allows controlling the strength of the connection between the sources X 1 and X 2 and the target Y .The simulation design showing the causal connections is reported in Fig. 1(a) alongside an example of the spectral content of each of the four processes when c = 1 (b).
The results of interaction information decomposition performed in the time and frequency domains for the VAR process (15) are shown in Fig. 2 (a, b-e, respectively).When the coupling parameter c is increased from 0 to 1, the coupling between the first source X 1 and the target Y increases, and the time-domain measure F X 1 ;Y (red line, Fig. 2 (a)) also increases as a direct consequence; an opposite trend is obtained for F X 2 ;Y (blue line, Fig. 2(a)) which decreases for increasing values of the coupling parameter c.The measure of source interaction (I X 1 ;X 2 ;Y , black line, Fig. 2(a)) is almost null when c = 0, while it highlights negative interaction information (denoting synergy) when c = 0.5, and finally becomes slightly positive (denoting redundancy) when c = 1.The frequency domain expansion of the information measures reveals interaction mechanisms which are specific for the oscillations simulated at 5 Hz and 35 Hz.Specifically, when c is low, there is a transmission of the oscillations from X 2 to the target through the link Z 3 → Y , with peaks at 5 Hz and 35 Hz in f X 2 ;Y and f X 1 X 2 ;Y (blue line, Fig. 2(b,d); the blue line in panel (b) is not visible as it is fully overlapped with the red line); the interaction information i X 1 ;X 2 ;Y is greater than 0 particularly at 35 Hz, denoting redundancy (Fig. 2(e)).When c = 1, there is a direct transmission of the oscillations at 5 Hz through the link X 1 → Y revealed by the peaks in f X 1 ;Y and f X 1 X 2 ;Y (red line, Fig. 2(b,c)); moreover, there is also a less prominent peak at 35 Hz due to the indirect causal link from Z 4 ∈ X 2 to the target, mediated by X 1 , which is characterized by a fully redundant interaction between the two sources (see the peak in i X 1 ;X 2 ;Y at 35 Hz, Fig. 2(e)).
When 0 < c < 1, the coexistence of redundancy, occurring at 35 Hz due to the cascade Z 4 → X 1 → Y , and synergy, occurring at 5 Hz due to the links X 1 → Y and Z 3 → Y , is made evident by analysing the spectral profile of the interaction information, but is not detectable employing the time-domain measure which indicates only a synergistic effect.In particular, when c = 0.5 (spectral profiles highlighted in green) the information shared individually between one source and the target (f X 1 ;Y , f X 2 ;Y ; Fig. 2(c,d)) shows comparable trends for the two sources, with a higher contribution at 35 Hz related with the oscillations of the common driver process Z 4 .On the other hand, the information shared jointly by the two sources and the target (f X 1 X 2 ;Y , Fig. 2(b)) still exhibits a higher contribution due to the autonomous oscillations at 5 Hz of Z 2 and Z 3 ; as a result, the interaction information (i X 1 ;X 2 ;Y , Fig. 2(e)) reveals greater synergistic contribution at 5 Hz if compared to the redundant contribution at 35 Hz.Importantly, the integrated measures in the time domain (Fig. 2(a)) suggest the impossibility of discriminating the coexistence of synergistic and redundant contribution with a theoretical value that is negative, indicating instead only a synergistic contribution (I X 1 ;X 2 ;Y < 0, Fig. 2(h)).

B. ESTIMATION
This section reports the practical estimation of the interaction measures defined in Sect.II.A,B, performed using the parametric and non-parametric approaches described in Sect.II.C, on simulated multivariate time series generated as realizations of the VAR(3) process (equations ( 15)).In this analysis we set the coupling parameter c = 0.5; this scenario, as previously seen with the theoretical example, guarantees the coexistence of redundant and synergistic information in two separate frequency bands.One hundred realizations of the processes were generated with a time-series length N = 1000, thus simulating a duration of 10 s with a sampling frequency of 100 Hz.
All the measures appearing in the interaction information decomposition were computed first identifying the VAR(3) model through the OLS method, followed by the computation of the Spectral Density matrix ( P Z , Eq. 6) from which all the time and frequency domain measures were extracted.With this approach, for each realization, the model order p was estimated through the minimum description length (MDL) criterion [39].The second approach consisted of the direct computation of all the terms of the Spectral Density matrix by using the WC method.Specifically, following [35], we choose the window suggested by Parzen, also known as the de la Vallèe-Poussin window, that shows a significantly lower side-lobe level compared to Hanning and Hamming windows; furthermore, it is non-negative for all frequencies, and produces non-negative spectral estimates [28].For the Parzen window, the relationship between the bandwidth (B w ) of the spectral window and the lag τ at which correlation estimates are truncated is B w = 1.273f s /τ .To resolve the corresponding peaks in the spectrum, we set the window bandwidth equal to 2.5 Hz, which brings to τ ≈ 51 using f s = 100 Hz.
To appraise the performance of parametric and nonparametric approaches, all the time and frequency domain measures included in the Interaction Information Decomposition of Eqs.(8,9) were computed as defined in Eqs.(7a,7b,7c,8,9) for each realization and method.Results of this computation are reported in Fig. 3, showing the distributions across the 100 realizations of all measures, compared with the exact values obtained from the true model parameters.The interaction information decomposition in the frequency domain highlights interaction mechanisms which are specific of the oscillations simulated at 5 Hz and 35 Hz, as in the previous subsection.The main spectral peaks displayed by all measures, whose theoretical profiles are shown by the red lines in Fig. 3, are well-resolved by both parametric and non-parametric estimation approaches.The bias and the variance of the spectral estimates are higher for the WC approach compared to a VAR model.This finding is confirmed when the measures are integrated over the whole frequency axis, as the non-parametric approach tends to overestimate all time domain measures, with a difference that is statistically significant for F X 1 ;Y , F X 2 ;Y and I X 1 ;X 2 ;Y (Wilcoxon test, p < 0.001).
To investigate the effects of an incorrect selection of the VAR model order in parametric estimation, the spectral profiles of the interaction information measure were computed for different values of the order in the range p ∈ {1, 2, 3, 10, 50}.In a similar way, the effects of windowing Comparison between parametric and non-parametric approaches for the estimation of time and frequency domain information measures computed over realizations of the four-node system (Eqs.15) simulated with coupling parameter c = 0.5.The figure depicts the distribution over 100 realizations (median and 5 th − 95 th percentiles) of spectral profiles (a-d) and time-domain values (e-h) of the information shared by the two sources and the target jointly (f X 1 X 2 ;Y , F X 1 X 2 ;Y (a, c) and individually (f X 1 ;Y , f X 2 ;Y , F X 1 ;Y , F X 2 ;Y , (b,c,f,g), and of the interaction information (i X 1 ;X 2 ;Y , I X 1 ;X 2 ;Y , (d,h)) estimated through parametric (PAR, blue lines) and non-parametric (No-PAR, orange lines) approaches.In each plot, the true values computed from the original model parameters are depicted as red lines (theoretical values for time domain analysis: F X 1 X 2 ;Y = 1.11,F X 1 ;Y = 0.35, F X 2 ;Y = 0.37, I X 1 ;X 2 ;Y = −0.38).Distributions are depicted as median and 5 th − 95 th percentiles of the frequency profiles, and as violin plots with probability densities estimated using the kernel density estimator (with dots representing outliers) for the time domain measures.Results of Wilcoxon test comparing PAR and No-PAR: *p < 0.005. in non-parametric estimation was investigated changing the spectral bandwidth in the range B w ∈ {30, 10, 5, 2.5, 0.5} Hz, which corresponds to truncating the correlation estimates at the lags τ ∈ {5, 15, 25, 55, 500}.The interaction information was computed for frequencies ranging from 0 to the Nyquist frequency (f s /2 = 50 Hz), and values indicative of the overall information shared among the processes in the low frequency (LF, 0.05 − 6.5 Hz), high frequency (HF, 34 − 36 Hz) bands, as well as over the entire frequency range (TOT, 0 − 50 Hz) were obtained.The LF and HF bands were centered around the two frequency peaks of the target process (5 Hz, 35 Hz), with bounds defined at the frequencies corresponding to −3 dB attenuation of the power spectrum.The average values of the interaction information within these bands, depicted as distribution of values over 100 realizations of the simulation, are depicted in Fig. 4. Looking at the estimates within the LF and HF bands (Fig. 4(a,b,d,e)), in the case of parametric estimation (blue violin plots) the interaction information exhibits the lowest bias when p = 3, which represents the true model order.The bias becomes very strong for model order p = 1, while it is much less evident when the model order close to the true value (p = 2) or is overestimated (p = 10, p = 50); the variance of the estimates tends to increase with higher model orders.In the case of the non-parametric approach (orange violin plots) the median value of the interaction information is quite far from the theoretical value when the number of lags considered for the estimation is low (τ = 5), and decreases progressively when the maximum lag used increases up to τ = 500; on the contrary, the variance of the estimates shows a clear tendency to increase with the truncation lag τ .
The distribution of the interaction information i X 1 ;X 2 ;Y (f ) integrated over the whole frequency axis, corresponding to the time domain measure I X 1 ;X 2 ;Y , is depicted in Fig. 4(c,f), and the corresponding estimates of bias and variance are reported in Table 1.The results confirm the opposite trends of the bias and the variance of the estimates as, for both estimation methods the bias increases when the variance decreases, and vice versa.

IV. APPLICATION TO PHYSICAL SCENARIOS
This section reports the application of the framework developed for the analysis of interaction between blocks of multivariate processes in the frequency domain.Effect of estimation parameters on spectral estimates of the interaction information of simulated time-series.Violin plots depict the distributions across 100 realizations of the interaction information I X 1 ;X 2 ;Y obtained integrating the frequency domain measure i X 1 ;X 2 ;Y (f ) computed using the parametric (blue box plot) and the non-parametric approach (orange box plot) over the low-frequency range (LF, 0.05-6.5 Hz; (a,d)), over the high-frequency range (HF 34-36 Hz, (b,e)) and over the whole frequency range (TOT, 0-50 Hz; (c,f)).The probability densities of the violin plots are estimated with the kernel density estimator; dots represent outliers, true values are indicated by the red dotted lines.

TABLE 1.
Values of bias and variance (average over 100 realizations of the simulation) of the time domain measure of interaction information obtained using parametric (PAR) and non-parametric (No-PAR) estimation approaches for different value of VAR model order p and truncation lag τ for the estimation of spectral and cross-spectral densities.
Three analyses were performed considering rather diverse complex systems which, perhaps counter-intuitively, present some notable common features.While brain units, electronic oscillators and climate systems exhibit oscillations in very different frequency ranges, a common thread throughout these systems is the fact that their subsystems are often synchronized across multiple time scales.In the brain, there are some dominant rhythms in different frequency bands which synchronize across spatial regions, often reflecting functional tasks such as memory, cognition and attention [40].Coupled electronic oscillators display common dynamics in time and space, which occur at different frequencies as a result of direct structural links but also non-trivial phenomena such as remote synchronization [41].Climate systems are influenced not only by land temperatures, but also by storage of heat and carbon dioxide in the oceans, which depends on physical and biological processes occurring across multiple time scales [42].
In the first application we analyzed the human brain, which provides a rich and important terrain to study the synchronization of neuronal populations and circuits endowed with higher-order interactions.In neuroscience, for example, a ubiquitous challenge consists of trying to infer functional brain networks from recorded EEG signals.However, this task, known as ''network inference'', is highly non-trivial and often reduced to the study of pairwise interactions [3].In this context, our first application concerns the analysis of high-order interactions as indexed by the EEG signals of 109 healthy subjects during the execution of a motor task [22].It is well-established that adults present 8-12 Hz (i.e.µ) and 16-26 Hz (i.e.β) rhythms in the EEG recorded from the scalp over the primary sensorimotor cortex.These rhythms show amplitude fluctuations that are synchronized to movement initiation or imagery [43].Besides specific oscillations in the frequency domain, a recent meta-analysis of functional MRI studies found contralateral primary motor cortex, contralateral somatosensory regions, and pre-motor regions activated during motor task execution [44].Given their clearly defined occurrence, the information derived from these oscillations and their localization on the scalp has also been used for developing brain-computer interfaces [45].
Subsequently, we considered a dataset of multiple time-series recorded from a unidirectionally-coupled ring wherein each unit consists of an electronic non-linear chaotic oscillator [46].These oscillators exhibit very rich dynamics that have previously been studied by considering pairwise interactions, estimating Granger causality and transfer entropy in their bivariate form [46], [47], however without considering higher-order interactions.Previous analyses have established a form of ''remote synchronization'' which takes place in this network of electronic oscillators, manifested as a dip and subsequent recovery of synchronization measures computed varying the distance between nodes in the ring [23], and have elucidated in detail the underlying mechanisms, consisting in amplitude modulation and interference phenomena between the higher and lower sidebands of a predominant oscillation in the system [46].The study of chaotic oscillators over the last years has raised considerable interest owing to the ability to replicate emergent properties characteristic of biological, particularly neural systems.For instance, the amplitude modulation and interference effects previously unveiled in the dataset analyzed here [46] have been previously observed in patterns of local field potentials of neuronal populations [48].Notably, an electronic circuit provides a system of drastically reduced scale and complexity compared to a physiological one, yielding full access to the activity of each node.This is akin to considering a micro-scale system and description, as opposed to a macro-scale one.
Finally, we applied interaction information decomposition in the frequency domain to discover higher-order interactions between time-series representative of the climate system.Knowingly, the Earth's climate system hosts a myriad of nested and coupled sub-systems, operating at diverse temporal and spatial scales [49].The climate system provided the initial inspirations for developing fundamental theories in chaos and complexity; moreover, it has been demonstrated that the climate dynamics are chaotic at specific scales while being stochastic or even regular at others [50].For example, different studies have proved the existence of cause-effect relationships between temperature and carbon dioxide concentration and how the latter can Granger-cause temperature increase [24], [38].For these reasons, higher-order measures in the frequency domain can be used to study the units of this system, showing how their oscillatory components can cooperate, when taken together or individually, in exchanging information.

A. EEG RECORDINGS
The analyzed dataset refers to EEG signals recorded, in 109 healthy participants, from 64 electrodes referenced to both mastoids as per the international 10-20 system with a sampling frequency of f s = 160 Hz [22].All data were drawn from a publicly-available dataset (https://physionet.org/content/eegmmidb/1.0.0/) and no new experiments were performed as a part of this work.The conditions analyzed in this work include two two-minute runs of a motor execution task performed as follows.First, a target was visualized on the right side of the screen (RIGHT condition), and the participants were asked to open and close the right fist cyclically until the target disappeared; afterwards, the participants relaxed (REST condition).The raw signals were firstly detrended, then filtered with two second-order Butterworth filters (band-pass, 2-35 Hz; notch, 59-61 Hz), and finally epoched to extract ∼23 trials for each condition and participant with a duration of 4 s each.For our analysis, we selected as sources X 1 and X 2 the groups of EEG channels [FC 2 , FC 4 , FC 6 ], [C 2 , C 4 , C 6 ] and as a target Y the group composed by [C 1 , C 3 , C 5 ] (Fig. 5 (a)).The data in each trial were tested for a restricted form of weak-sense stationarity [51] and, subsequently, parametric estimates of the spectral information measures were obtained for each subject, experimental condition and trial.A VAR model was identified from the nine selected time-series setting the model order according to the MDL criterion.The Durbin-Watson test for whiteness of the innovations, testing the absence of serial correlations in the VAR residuals, was performed [33]; then, the estimated VAR parameters (VAR coefficients and input covariance matrix) were used to compute the PSD matrix for the three multivariate processes according to (6), and the frequency-domain functions measuring the information shared between the target and the two sources and the interaction information were obtained as described in Section II.Moreover, to assess the statistical significance of interaction information in the frequency domain, the surrogate data analysis described in Sect.II-D was applied separately for each trial, participant and condition.Finally, values indicative of the overall information shared within the α and β bands were obtained by integrating the measures over the appropriate frequency ranges; e.g., for the interaction information, I X 1 ;X 5(b,c,d,e) reports the grand-average over participants and trials of the frequency profiles for each information measure, computed separately for the REST and RIGHT conditions and depicted over the frequency range under consideration (1-40 Hz).For each frequency-domain measure, the spectral profiles feature a prominent peak around 10 Hz (α band), almost unchanged between the two conditions, and a reduction of the average value in the β band (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26) during the RIGHT compared to the REST condition.
The spectral measures quantified statistically significant amounts of information shared between the processes at each frequency.In particular, in Fig. 5(f,g) the grand-average profiles of i X 1 ;X 2 ;Y (f ) computed in the two conditions are compared with the 97.5 th percentile of the distribution of the measure computed for the surrogate time series, showing a statistically significant redundant contribution in the network with a prominent peak in the α band, which tends to decrease towards 30 Hz. (b-e) Spectral profiles (mean (bold line) and 25 th − 75 th percentiles (shades) over 109 participants) of the information shared between the two sources and the target jointly (f X 1 X 2 ;Y ) and individually (f X 1 ;Y , f X 2 ;Y ), and of the interaction information (i X 1 ;X 2 ;Y ) obtained during relaxation (REST, green lines) and while closing the right fist (RIGHT, orange lines).(f-g) Distribution of the spectral profiles of the interaction information computed separately for REST and RIGHT conditions, compared with surrogate distributions (gray lines).(h-i-l) Violin plots depicting the distribution across participants of the information measures obtained by integrating the corresponding spectral measures over the whole spectrum (h), within the α band (8-12 Hz, (i)) and within β band (16-26 Hz, (l)), computed during relaxation (REST) and closing the right fist (RIGHT).Wilcoxon test comparing REST and RIGHT conditions: * , p < 0.05, corrected for multiple comparison across four measures.
Fig. 5(h,i,l) depicts the violin plots for the distributions across participants of all spectral information measures integrated over the whole frequency spectrum, α and β bands, computed separately for the REST and RIGHT conditions.By analyzing the distributions obtained for each measure, the information shared by the two sources taken together and the target (F X 1 X 2 ;Y ) assumes the highest values, generally between 1 and 3 nats, the amounts of information shared by the target and a single source are comparable (with a slight prevalence of F X 2 ;Y over F X 1 ;Y ), and the lower values are attained by the interaction information I X 1 ;X 2 ;Y .A tendency towards a decrease of all information measures is noticeable during motor task execution (RIGHT, orange violin plots) compared to the resting condition (REST, green violin plots).Such a decrease is not significant for any measure in the whole spectrum and in the α band, while it is statistically significant (Wilcoxon test, p < 0.05) for the measures F X 1 ;Y and F X 1 X 2 ;Y computed in the β band.
Thus, the application of our framework to real EEG signals recorded during the motor execution task showed a statistically significant reduction of the information exchanged in the β frequency range (16-26 Hz) associated with the execution of the movement of the right hand.This result can be related with the physiological phenomenon of event-related desynchronization [52].Importantly, the time-domain measures and the measures computed in the frequency range where the EEG activity is prominent (α band, 8-12 Hz) did not vary significantly during the task.These findings highlight the importance of analyzing multivariate EEG interactions within specific frequency bands to obtain markers of motor execution potentially useful for the design of brain-computer interfaces.

B. ELECTRONIC CHAOTIC OSCILLATORS
The second application is relevant to a ring of coupled electronic oscillators, for which the structural diagram of the oscillator circuit is visible in Fig. 6(a  sampled with a sampling frequency f s = 100 kHz and are freely available [53]. The frequency spectrum of each node presents three distinct peaks: the most prominent (central one) at f c ≈ 2.8 kHz, and two weaker ones (sidebands) at f l = f c /2 ≈ 1.4 kHz and f h = f l + f c ≈ 4.2 kHz.The higher sideband represents the mirror frequency of the lower one.Due to the presence of the integrators and their associated saturation non-linearity, another oscillatory component is generated in the frequency range of the lower sideband, caused by the demodulation action of the system.Depending on the phase relationship between the low-frequency component of the input signal and such demodulated signal, constructive or destructive interference can arise [46].These phenomena, taken together, lead to spatial fluctuations of the lower sideband amplitude, which engender a pattern with features closely related to the remote synchronization effect [41].In this system, remote synchronization manifests as a non-monotonic decay of synchronization along the ring, wherein the average synchronization drops with increasing distance from a given node, then increases transitorily, and finally vanishes [46].
As said, it is known that in this system, moving away from a node, synchronization initially decays, then gradually increases up until a distance ≈8 nodes, and eventually vanishes [46], [47].However, the structural couplings on the ring are only between first neighbors as indicated by the master-slave configuration in Fig. 6(a).For these reasons, we choose the target and sources to be analyzed with multivariate information measures as reported in Fig. 6(c); specifically, denoting {S i } 32 i=1 the set of voltage signals recorded from each oscillator in the ring, the target process is varied alongside the ring structure (Y = S i ), the first source is kept at a distance of two nodes from the target (X 1 = S i+2 ), and the second source is varied keeping a distance from the target of at least four nodes (X 2 = S i±d , d ≥ 4).Before the analysis, for each network node i, the three time-series measured as voltage output by the two sources and the target were sub-sampled by a factor of 10 (yielding a sampling frequency of f s = 10 kHz) and normalized to zero mean.Then, for each node, the spectral measures and their integrated versions in the time domain were obtained as follows.First, the data were tested for a restricted form of weak-sense stationarity [51]; Then, a VAR model was identified on the three time-series using the OLS method to estimate the model coefficients and setting the model order according to the MDL criterion.Afterwards, the Durbin-Watson test for whiteness of the VAR residuals was performed [33], the estimated model parameters were used to compute the PSD matrix and the information measures for frequencies ranging from 0 to the Nyquist frequency f s /2 = kHz, the overall shared among the processes was determined integrating the measures over entire frequency range.Furthermore, to assess the statistical significance of interaction information in the frequency domain, the analysis of surrogates as described in II-D was performed.The entire analysis was then repeated considering, for any given target node S i , i = 1, . . ., 32, each of the 25 possible triplets {S i , S i+2 , S i±d }, d ≥ 4.
The results of the interaction information decomposition in the frequency domain are reported in Fig. 7(a-d).The profile of each measure is evaluated for each target and each distance d, and then averaged over all targets in the ring.In general, all analyzed measures reveal information shared among the processes at frequencies corresponding to the central peak (f = 2.8 kHz) and to the two side-bands (lower, f = 1.4 kHz and higher, f = 4.2 kHz).The information shared between the source X 1 and the target Y (distance from the target equal to 2) reported in Fig. 7(a) shows completely overlapped trends regardless of the distance between X 2 and Y .This is not the case for the information shared between X 2 and Y and the information shared jointly between {X 1 , X 2 } and Y , which show a reduction of the corresponding information measures as the distance between X 2 and the target increases (Fig. 7(b,c)); this is particularly evident in the frequency ranges corresponding to the lower and higher side-bands.The analysis of interaction information reported in Fig. 7(d) reveals a strong redundant contribution at ≈3 kHz corresponding to the central peak, as well as  (a-d) Spectral profiles of the information measures (f X 1 ;Y , f X 2 ;Y , f X 1 X 2 ;Y and i X 1 ;X 2 ;Y ) averaged over the 32 target processes for the entire ring.The profiles are reported for four different distances of the source X 2 from the target, with d = 4 (blue line), d = 8 (black line), d = 12 (green line) and d = 16 (red line).(e) Spectral profile of the interaction information for a distance d equal to 4 and averaged along with the 32 ring nodes (blue line), and corresponding 97.5 th percentile computed across 100 surrogates averaged over the 32 oscillators in the ring (orange line).(f) Spectral profiles of interaction information integrated over the entire frequency range for different distances between the considered target i and the source X 2 .
a synergistic contribution of both the lower and higher side-bands for a distance equal to 4, which gradually vanishes with the increasing of the distance (Fig. 7(d)).All these interactions, including the synergistic ones found in the two sidebands for d = 4, resulted as statistically significant according to the surrogate data analysis (Fig. 7(e)).
The analysis of the integrated measures (Fig. 7(f)) highlights that the information shared between the source X 1 and the target Y is independent of the distance d between the target and the second source X 2 (red line), while both the information shared between X 2 and Y (blue line) and the information shared between {X 1 , X 2 } and Y (green line) decrease with the distance, reaching a minimum for d ≈ 16 (farthest point from the target).Moreover, the interaction information (black line) shows the prevalence of redundancy between the three processes for all distances, with a maximum value at d ≈ 8 and a minimum at d ≈ 16.
In summary, this section sought clarity on the mechanisms underlying the constructive or destructive interference in the ring that is mostly related with the interaction between lower and higher sidebands.This has been performed in both time and frequency domains by analyzing the presence of redundancy or synergy in the system.By doing so, we found that the highest value of redundancy is related to the most prominent peak of activity (f = 2.8 kHz), but also that synergy can arise depending on the distance between X 2 and the target.In particular, significant synergistic effects emerge for d = 4, reflecting coupling mechanisms between oscillators not directly connected via physical links.These findings elucidate the complex nature of the interactions between the analyzed oscillators, and suggest that the linear analysis of high-order interactions confined to frequency bands related to known demodulation effects [46] can disclose constructive synchronization phenomena in nonlinear systems.

C. CLIMATE DYNAMICS
In the third application we considered the interaction information decomposition of time-series representing the dynamic activity of three environmental variables: global land temperature (TL), global temperature of the ocean (TO), and carbon dioxide volume (CO 2 ).These three variables are strictly interconnected with the carbon cycle.In particular, the latter is consumed mainly through photosynthesis on land and by absorption in cold ocean waters, while it is released from warm ocean waters; furthermore, the circulation of ocean waters from warm to cold zones, and vice versa, promotes both consumption and production of CO 2 , with rates that are strictly dependent on temperature and gas concentration [54].The TL and TO time-series used in this study were downloaded from the NOAA National Center for Environmental Information [55], providing global temperatures in terms of the average land and ocean temperature measured over the 20 th century.The CO 2 time-series was collected by the Mauna Loa Observatory, Hawaii [56] and consists of the monthly average mole fraction of CO 2 in the atmosphere, i.e., the number of molecules of CO 2 in every one million molecules of dry air.Dates starting from January 1980 to April 2019 (maximum overlap between the CO 2 and temperature datasets) were selected, so that the resulting dataset includes 3 time-series comprising 492 points each.Each time-series was first detrended by applying an l 1 -norm filter to fulfill stationarity requirements, and then normalized to zero mean and unit variance.Subsequently, the seasonal oscillatory component was removed by calculating the individual monthly temperatures and subtracting them from the time-course temperatures for that month.The resulting time series are shown in Fig. 8(a-c).The data were then tested for a restricted form of weak-sense stationarity [51] In this case, we performed non-parametric estimation of the measures of information decomposition computed in both time and frequency domains.Specifically, given that data have a sampling period of 1 month, we selected a Parzen window with bandwidth B w = 0.025 month −1 , which corresponds to a number of truncation points of the correlation functions τ = 51 (see Eq.13).The information measures were derived selecting the global land temperature as target process (Y =TL), and the temperature of the ocean and the carbon dioxide volume as sources (X 1 = TO X 2 = CO 2 ).The statistical significance of the interaction information measure computed in the frequency domain was performed via the method described in II-D, setting the same τ for each of the 200 surrogate triplets.
Fig. 8 shows the spectral profiles of information shared by the TL (target, Y ) and the two sources (X 1 = TO, X 2 = CO 2 ) (Fig. 8(d)) and of the interaction information among the three processes (Fig. 8(e)); the latter is depicted together which its distribution over 200 surrogate triplets.The interaction information show a significant redundant contribution (i X 1 ;X 2 ;Y > 0) at low frequencies (corresponding to a time scale of about 32 months) that is mostly due to the similar trends of f X 1 ;Y (red line) and f X 1 X 2 ;Y (green line).On the other hand, the measure computed at high frequencies (time scale of about 4 months) reveals a statistically significant synergistic contribution (i X 1 ;X 2 ;Y < 0), indicating that the two sources provide additional information about the target when they are considered jointly rather than separately.The time domain analysis returns instead the following integrated measures: F X 1 ;Y = 0.13, F X 2 ;Y = 0.06, F X 1 X 2 ;Y = 0.16 and I X 1 ;X 2 ;Y = 0.02.Interestingly, the integrated interaction measure (I X 1 ;X 2 ;Y ) takes a small positive value suggesting the presence of net redundancy, thus preventing the detection of coexisting redundancy and synergy which is instead revealed at different time scales looking at the spectral measure.
In summary, the results of our analysis of climate data highlight again the importance of measuring higher-order interactions in the frequency domain.Indeed, while the time-domain analysis of high-order interactions indicates only a small amount of redundancy, and while pairwise interactions in the frequency domain merely reflect the coupling between TL and TO and between TL and CO 2 at ∼4 and ∼32 months, the spectral measure of interaction information reveals non-trivial mechanisms at the shorter time scales.The main property disclosed by this measure is the emergence of synergy at ∼32 months, suggesting that the land temperature , pre-processed for information analysis.(d) Spectral profiles of the information shared between the target (Y = TL) and the two sources (X 1 = TO, X 2 = CO 2 ) taken jointly (f X 1 ,X 2 ;Y , green line) and individually (f X 1 ;Y and f X 2 ;Y , red and blue lines).(e) Spectral profile of the interaction information (i X 1 ;X 2 ;Y , blue line), plotted together with the distribution of the measure computed over 200 surrogate time series (2.5 th , 50 th and 97.5 th percentiles, brown line and green shaded area).
is linked to ocean temperature and carbon dioxide volume through different mechanisms at this time scale.

V. DISCUSSION
The present work extends the framework recently proposed for the frequency-domain evaluation of high-order interactions among stochastic processes [19] to the analysis of network systems formed by an arbitrarily large number of nodes.Moreover, the framework is endowed with a model-free estimation approach that is compared with the parametric estimator in simulated data.Finally, it is validated experimentally in three real-world scenarios relevant to network systems that produce output signals oscillating in very different frequency ranges and displaying signatures of pairwise coupling as well as higher-order effects.

A. THEORETICAL EXAMPLE AND SIMULATION
The theoretical example was designed to show how redundant and synergistic interactions between oscillations manifested with the same frequency in different signals can arise from the interplay between direct and indirect coupling mechanisms, and can be elicited from the spectral profiles of the information measures better than from their time domain values.We showed that redundancy may arise as a consequence of indirect interaction effects, such as the coupling between a source and the target mediated by another source, synergy may arise from interactions occurring independently between each source and the target at nearby frequencies, and that synergy and redundancy may coexist being simultaneously present in different frequency bands (Fig. 2(e)).Importantly, these frequency-specific interaction mechanisms may be hidden if one looks at the standard time-domain representation of information, as a result of the integration between positive and negative interaction information values produced at different frequencies (Fig. 2(e) vs. Fig.2(a)).The results obtained here in terms of the time-domain interaction information are in agreement with the observation of synergy and redundancy at different time scales made in Ref. [38] by employing a more sophisticated approach based on state-space models and partial information decomposition [14].
The following simulation study was conducted to assess the performance of parametric and non-parametric estimation approaches.First, we focused on quantifying the bias and variance of the two estimators in both time and frequency domains.The results shown in Fig. 3 highlight better performance, in terms of lower bias and variance, of the parametric approach compared with the non-parametric one.This is not surprising as the time series were generated from a VAR process with known model order; at any rate, the non-parametric estimator displayed a good capability to replicate the theoretical profiles of the information measures.The second analysis aimed to show the dependence of the interaction information measure on the estimation parameters.The results reported in Fig. 4 indicate that both the parametric and the non-parametric approach exhibit considerable variability of the interaction information with the estimation parameters (respectively, the model order p and the length τ of the window used to compute correlation sequences).Both approaches are subject to the bias-variance trade-off [31]: in parametric estimation, underestimation and overestimation of the order p lead respectively to increased bias and increased variance; in non-parametric estimation large values of τ allow the spectral window to behave as an impulse function, leading to smaller bias and higher variance.Furthermore, a non-parametric approach requires proper selection of the window type and bandwidth (inversely related to τ ) in order to obtain the desired trade-off between accuracy and spectral resolution.Both methodologies have merits and drawbacks: non-parametric estimates are more generally applicable, as the WC estimator does not make assumptions about the data except broad-sense stationarity; parametric estimates, being based on a rational transfer function model of the data, are more restrictive in this sense, but yield more accurate spectral estimates (and thus information estimates) when the model fits appropriately the data.On the other hand, when the model is incorrect, non-parametric approaches outperform the parametric ones [31].As a practical hint, we recommend testing the model assumptions (i.e., whiteness and uncorrelation of the model residuals) after model identification, and evaluate the utilization of classical spectral estimates when such assumptions are not verified.

B. EEG RECORDINGS
The study of human brain activity during motor imagery or motor execution plays an important role in clinical neurology and neuroscience.Any motor action results from a dynamic interplay between separate brain regions and involves diverse processes underlying movement preparation and execution.A well-established body of literature indicates an intrinsic balance between excitatory and inhibitory couplings among brain regions and between the hemispheres [57].According to this knowledge, in our multivariate analysis we selected the three EEG signals recorded from channels located on the hemisphere contralateral to the right-hand motor execution (Y = [C 1 , C 3 , C 5 ]) as vector target process, and two different groups of EEGs from the ipsilateral hemisphere as sources ( With this choice, the activity of the target could be related to the primary motor cortex (M1) of the contralateral (left) hemisphere, and that of the sources to the prefrontal motor cortex (PMC) and the supplementary motor area (SMA) of the ipsilateral (right) hemisphere [58].
Functional MRI studies have demonstrated that M1, SMA, and PMC are strongly involved in the planning and execution of hand motor tasks [59], and identified a well-defined brain network supporting the execution of right-hand motor execution tasks [60], [61].Specifically, a positive coupling parameter has been inferred between M1, SMA, and PMC of the contralateral hemisphere, accompanied by cross-inhibition of neural activity between the two M1 areas during task execution.Our results show a reduction of the measure of information shared between the EEG activity recorded from the electrodes associated with M1 and PMC (F X 1 ;Y ), as well as with M1 and {PMC,SMA} (F X 1 X 2 ;Y ), during motor execution compared to the REST condition (Fig. 5).Notably, the reduction turns out to be statistically significant only when the information measures were integrated within the β frequency range (16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26).This modification of the information shared between groups of EEG signals could be related to event-related desynchronization effects [52].
Our findings are supported by other studies performing brain connectivity estimation in both time and frequency domains, which documented a decreasing degree of connectivity associated with motor task execution in the contralateral hemisphere [8], [62], [63].Nevertheless, while in this application we demonstrated the usefulness of evaluating information measures within specific frequency bands, we did not find distinctive patterns of high-order interactions showing the potential to improve the performance of BCI systems beyond the use of classical event-related measures [52].The dominant redundancy detected in all experimental conditions is likely related to effects of volume conduction that blur connectivity patterns identified at the level of scalp EEG sensors [64].For more informative analyses also allowing a more precise biological interpretation of the inferred connectivity patterns, we recommend the computation of spectral information measures at the level of the cortical sources identified applying source reconstruction techniques [65].

C. ELECTRONIC CHAOTIC OSCILLATORS
The time-series considered in the second application index the dynamic activity of a ring of thirty-two chaotic oscillators, for which the master-slave unidirectional structure guarantees a higher level of stationarity and more elementary dynamics with an a-priori well-known topological structure compared to physiological systems.The literature pointed out how systems like this replicate the spontaneous formation of multi-scale community structure as a function of coupling strength, with similarities to the modular organization observed of brain networks [66], [67].For this reason, it is reasonable to assume that chaotic oscillators could be taken as a benchmark for testing in physical scenarios new methods developed to study the interactions between more intricate dynamical systems.
In this context, the application was devised to test the proposed framework in a system where the emergence of complex behaviors was previously rigorously explained, and could support the interpretation of the patterns of multivariate information assessed in the frequency domain.In particular, in this specific system the non-linearity of the oscillators enables substantial energy exchange between the lower and higher side-bands (1.4 kHz and 4.2 kHz) through demodulation and interference phenomena [46].The initial experimental investigations of this system [67] have identified that the point of maximally-destructive interference between the lower side-band and the demodulated signal occurs at a distance between nodes around six steps [46].These complex nonlinear interaction and demodulation effects manifested as an interplay among the three different peaks are reflected in our frequency specific analysis of interaction information by the co-occurrence of redundant and synergistic information exchange.According to our results (Fig. 7), redundancy is naturally related to the information carried by the main oscillation at 2.8 kHz, while synergy is revealed for the lower side-bands mirroring the interference phenomena observed in [46].Also, the distance between nodes is in agreement with the detection of remote synchronization, as the statistically significant synergistic behaviors are associated with the lower side-band when d < 6 and become entirely redundant for d > 6. Contrary to the insights provided by frequency-specific analysis, the time domain measure obtained integrating the interaction information across all frequencies depicts purely redundant contributions (Fig. 7(f)).This finding, though being far less informative about the coupling mechanisms arising in the system, is in agreement with previous investigations [41], [66].In fact, the peak of the interaction information detected at a distance d = 8 between the second source and the target is consistent with a maximum of the cross-correlation coefficient at a distance d ≈ 8 over the ring between node pairs not directly connected with a physical link but connected with a group of intermediary nodes weakly synchronized with them [66].

D. CLIMATE DYNAMICS
As a third application, we studied the interactions between three time-series representative of the behavior of the global climate system: the global temperature of the land, the global ocean temperature and the carbon dioxide concentration.These time-series have been widely studied using bivariate analysis, documenting a strong correlation between global sea surface temperature variability and land surface temperature [68] as well as between the rise of CO 2 and global warming [69].Previous analyses based on linear modelling have focused on the inference of directed links, demonstrating that the level of CO 2 in the atmosphere causes the global land temperature [24], [70], but also showing that higher global temperatures promote a rise of greenhouse gas levels, thus establishing a positive feedback [71]; the direction of interaction between land temperature and CO 2 appears to be dependent on the time-scale [38].
The observation of coupling patterns dependent on the time scale motivates our frequency-specific analysis of interaction effects.Specifically, we found that the circuit formed by global land and ocean temperatures and CO 2 emissions is significantly redundant at low frequencies (time scale ∼32 months), and significantly synergistic at higher frequencies (time scale ∼4 months).The redundancy measured at the longer time scales can be related with the positive feedback between CO 2 concentration and the average temperature [54], documenting common mechanisms underlying the variations of land or ocean temperature and the dynamics of CO 2 .On the other hand, the synergistic effects observed at the faster scale (∼4 months) can be related with the turning of the seasons, which seems to involve separate mechanisms in driving the coupling between land and ocean temperatures on the one side, and land temperatures and CO 2 on the other side.The fact that the synergistic mode is found to be faster than the redundant mode, in this example, is also worth of further investigation to check whether it corresponds to a general feature of these systems and if the same phenomenon holds in other data sets.

VI. CONCLUSION AND IMPLICATIONS
Collective phenomena in complex network systems often emerge from multiple links among many elementary subsystems that may occur through different mechanisms and can be detected only going beyond the framework of VOLUME 9, 2021 pairwise interactions.This work demonstrates that higherorder interactions in a multivariate set of dynamic processes can remain hidden if they are investigated exclusively in the time domain, and that such processes may display redundancy in a certain frequency range and synergy in another range.This implies that in some applications of network measures to real-world systems the frequency decomposition of information-theoretic measures is mandatory.We documented this aspect in three rather diverse systems endowed with rich oscillatory dynamics, obtaining results that can foster the implementation of practical applications in the relevant field of study.
In the study of collective brain dynamics where information is typically retrieved within specific frequency bands, we showed the importance of focusing on the β EEG waves while analyzing multivariate interactions during motor execution and imagery tasks.While our results provide a discriminative parameter that can serve as an additional feature to be employed in brain-computer interfaces, the identification of non-pairwise interactions at the level of individual brain rhythms may also aid in the recognition of pathologically altered brain states such as schizophrenia, epilepsy, or dementia [38], [72].
In the study of coupled chaotic oscillators, the inference of pairwise and high-order interactions typically resorts to nonlinear methods [73].However, the results relevant to the ring of chaotic electronic oscillators analyzed in our second application showed that the effect of remote synchronization [23] can be described also using linear spectral analysis, which is more easily applicable on short and noisy experimental time series.This opens up new possibilities for extracting measurements from collective activity in distributed sensing applications.
Finally, in the analysis of climate dynamics, the classical modeling approaches [74], [75] are nowadays complemented by approaches relying solely on observations [76], which provide evidence for a statistically meaningful relation among anthropogenic emissions, natural cycles, and global temperatures.
To our knowledge, this is the first time that higher order statistical effects are analyzed in the coupled system of temperature and greenhouse gases; we showed that a given set of variables derived from the climate system may display redundancy and synergy over widely separate frequency ranges.Though preliminary, these results may bring additional knowledge useful to elucidate the mechanisms related to the large-scale drivers of regional climate change, and may possibly contribute to categorize the effects of human activity on global climate and, in perspective, to plan strategies to counteract such effects.
where | • | stands for matrix determinant.The quantity in (3) is a symmetric, non-normalized measure of coupling, which takes its minimum value f X ;Y (ω) = 0 when X and Y are uncorrelated at the frequency ω (|P [XY ] (ω)| = |P X (ω)||P Y (ω)|), and tends to infinity when X and Y are fully dependent at the frequency ω (|P [XY ] (ω)| = 0).Integration of (3) over the whole frequency axis yields the measure p n = [Z n−1 , . . ., Z n−p ] T and considering N consecutive time steps, a compact representation of the VAR model (11) can be defined as z = Az p + U , where A = [A 1 , . . ., A p ] is the M × pM matrix of unknown coefficients, z = [Z p+1 , . . ., Z N ] and U = [U p+1 , . . ., U N ] are the M × (N − p) matrices and z p = [Z p p+1 , . . ., Z p N ] is a pM × (N − p) matrix collecting the regressor terms.

FIGURE 1 .
FIGURE 1. (a)Graphical representation of the four-variate VAR process generated according to(15), where the network nodes represent the simulated scalar processes and the arrows represent the imposed causal interactions (self-loops depict the influences from the past to the present sample of a process).(b) Power spectral density of the four processes obtained with c = 1, revealing a prominent activity of both sources and target at 5 Hz and 35 Hz.

FIGURE 2 .
FIGURE 2. Time-and frequency-domain multivariate information measures computed for the simulated VAR process of Fig.1.(a) time-domain measure of information shared jointly and individually between the two sources and the target (F X 1 X 2 ;Y and F X 1 ;Y , F X 2 ;Y ) and of the interaction information (I X 1 ;X 2 ;Y ) plotted as a function of the coupling parameter c. (b-e) spectral profiles for the four measures plotted as a function of the frequency (f (Hz))for different values of coupling parameter (c = 0, blue lines; c = 0.5, green lines; c = 1, red lines).

FIGURE 3 .
FIGURE 3.Comparison between parametric and non-parametric approaches for the estimation of time and frequency domain information measures computed over realizations of the four-node system (Eqs.15) simulated with coupling parameter c = 0.5.The figure depicts the distribution over 100 realizations (median and 5 th − 95 th percentiles) of spectral profiles (a-d) and time-domain values (e-h) of the information shared by the two sources and the target jointly (f X 1 X 2 ;Y , F X 1 X 2 ;Y (a, c) and individually (f X 1 ;Y , f X 2 ;Y , F X 1 ;Y , F X 2 ;Y , (b,c,f,g), and of the interaction information (i X 1 ;X 2 ;Y ,

FIGURE 4 .
FIGURE 4. Effect of estimation parameters on spectral estimates of the interaction information of simulated time-series.Violin plots depict the distributions across 100 realizations of the interaction information I X 1 ;X 2 ;Y obtained integrating the frequency domain measure i X 1 ;X 2 ;Y (f ) computed using the parametric (blue box plot) and the non-parametric approach (orange box plot) over the low-frequency range (LF, 0.05-6.5 Hz; (a,d)), over the high-frequency range (HF 34-36 Hz, (b,e)) and over the whole frequency range (TOT, 0-50 Hz; (c,f)).The probability densities of the violin plots are estimated with the kernel density estimator; dots represent outliers, true values are indicated by the red dotted lines.

FIGURE 5 .
FIGURE 5. (a) Overview of the EEG electrode montage highlighting the positions of the target process Y covering the he contralateral motor area (channels [C 1 , C 3 , C 5 ]) and of the two source processes X 1 and X 2 , covering the ipsilateral premotor areas (channels [FC 2 , FC 4 , FC 6 ] and ([C 2 , C 4 , C 6 ]).(b-e) Spectral profiles (mean (bold line) and 25 th − 75 th percentiles (shades) over 109 participants) of the information shared between the two sources and the target jointly (f X 1 X 2 ;Y ) and individually (f X 1 ;Y , f X 2 ;Y ), and of the interaction information (i X 1 ;X 2 ;Y ) obtained during relaxation (REST, green lines) ).The circuit consists of four summing stages associated with low-pass filters.Three of such stages with negative gains G 1 = −3.6,G 2 = −3.12,G 4 = −3.08 and filter frequency F 1 = F 2 = F 3 = 2 kHz are arranged as a ring oscillator.Two Integrator stages with integration constants K 1 = 3.67, K 2 = 0.11µs −1 and mixing gains G 3 = −0.5 and G 5 = −0.71are overlapped to this structure.The ring is completed by the fourth summing stage having F 4 = 100 kHz F 1 with one input (gain G 6 = 0.132) which is necessary to close the internal ring itself and another input (gain G i = −1.44)connected to the previous oscillator in the ring network (Fig.6(b)).To limit the voltage swing, an inverter with gain G 0 = −0.4 is provisioned.The recorded time-series have a length l = 65536 points, are

FIGURE 6 .
FIGURE 6.(a) Diagram of the oscillator circuit instanced at each network node.(b) Master-slave (unidirectional, clock-wise) structure of the ring comprising 32 oscillators.(c) Graphical representation of the three-variate system analyzed for each target system S i , i = 1, . . ., 32.

FIGURE 7 .
FIGURE 7.(a-d) Spectral profiles of the information measures (f X 1 ;Y , f X 2 ;Y , f X 1 X 2 ;Y and i X 1 ;X 2 ;Y ) averaged over the 32 target processes for the entire

FIGURE 8 .
FIGURE 8. (a-c) Normalized time-series of global land Temperature (TL), global temperature of the ocean (TO), and carbon dioxide volume (CO 2 ), pre-processed for information analysis.(d) Spectral profiles of the information shared between the target (Y = TL) and the two sources (X 1 = TO, X 2 = CO 2 ) taken jointly (f X 1 ,X 2 ;Y , green line) and individually