# A whitening approach for Transfer Entropy permits the application to narrow-
band signals

Daube, Christoph1∗
Gross, Joachim1,2
Ince, Robin A. A.1

( 1Institute of Neuroscience and Psychology, University of Glasgow, 62
Hillhead Street, Glasgow G12 8QB, UK
2Institute for Biomagnetism and Biosignalanalysis, University of Münster,
Malmedyweg 15, 48149 Münster, Germany
*Correspondence<EMAIL_ADDRESS><EMAIL_ADDRESS>)

## Abstract

Transfer Entropy, a generalisation of Granger Causality, promises to measure
“information transfer” from a source to a target signal by ignoring self-
predictability of a target signal when quantifying the source-target
relationship. A simple example for signals with such self-predictability are
narrowband signals. These are both thought to be intrinsically generated by
the brain as well as commonly dealt with in analyses of brain signals, where
band-pass filters are used to separate responses from noise. However, the use
of Transfer Entropy is usually discouraged in such cases. We simulate
simplistic examples where we confirm the failure of classic implementations of
Transfer Entropy when applied to narrow-band signals, as made evident by a
flawed recovery of effect sizes and interaction delays. We propose an
alternative approach based on a whitening of the input signals before
computing a bivariate measure of directional time-lagged dependency. This
approach solves the problems found in the simple simulated systems. Finally,
we explore the behaviour of our measure when applied to delta and theta
response components in Magnetoencephalography (MEG) responses to continuous
speech. The small effects that our measure attributes to a directed
interaction from the stimulus to the neuronal responses are stronger in the
theta than in the delta band. This suggests that the delta band reflects a
more predictive coupling, while the theta band is stronger involved in bottom-
up, reactive processing. Taken together, we hope to increase the interest in
directed perspectives on frequency-specific dependencies.

## Introduction

Over the last decades, the description of statistical dependencies in cerebro-
cerebral and cerebro-peripheral pairs of time series has witnessed a surge of
interest (Bassett & Bullmore, 2006; Brookes et al., 2011; Naselaris et al.,
2011; Crosse et al., 2016; Mell et al., 2021; Gross et al., 2021). In these
fields, the general idea is to gain insight into the workings of the brain by
either studying how time series of neuronal activity relate to other neuronal
activity or to external signals such as auditory or visual stimuli as well as
the activity of other organs.

As a consequence, countless methodological approaches have been suggested to
mathematically quantify these dependencies (Bastos & Schoffelen, 2016). Some
of these ideas specifically aim at the description of directed interactions,
for example by using measures of the so-called “Granger-causal” (Granger,
1969) family, or their generalisation to nonlinear relationships, Transfer
Entropy (Schreiber, 2000; Barnett et al., 2009). In both of these measures,
the main idea is to quantify a directional dependency by first assessing to
what degree a target time-series can be predicted from itself and secondly
assessing to what degree this auto-prediction can be improved upon with the
assumed source time-series. In its classic formulation, TE implements this by
way of conditioning the mutual information (MI) between source and target on
an operationalisation of the target past. It has been suggested that this
warrants the capacity to correctly estimate not only “predictive information
transfer” between source and target time-series, but also the recovery of the
true underlying interaction delay (Wibral et al., 2013). The precise
estimation of such quantities is of high interest for the research programmes
not only of functional connectivity, but of cognitive neuroscience in general.

However, in the arguably simplest and in many applications ubiquitous case of
self-predictable or auto-correlated time-series, namely narrowband time-series
as obtained e.g. when applying band-pass filters, TE fails to deliver
intuitively comprehensible results. The application of TE in such cases has
therefore repeatedly been discouraged (Florin et al., 2010; Barnett & Seth,
2011).

Such applications however are of potentially high interest, given that
frequency specific interactions are at the core of popular hypotheses about
cerebro-peripheral (Giraud & Poeppel, 2012; Donhauser & Baillet, 2020) and
cerebro-cerebral (Schnitzler & Gross, 2005; Fries, 2015; Michalareas et al.,
2016; Schoffelen et al., 2017) interactions. There is an abundance of evidence
for intrinsically auto-correlated or band-limited parts of neuronal activity
(“oscillations”), whose presence in neuronal recordings (Wang, 2010; Donoghue
et al., 2020) should accordingly impede the use of TE even without the use of
analysis filters. Moreover, in light of the usually low signal-to-noise ratio
(SNR) of many recording modalities, particularly of non-invasive neuroimaging,
the isolation of band-limited activity via spectral filtering is a pervasive
strategy to achieve acceptable sensitivity and specificity. The main idea of
TE thus turns out to be an empty promise for many real-world applications
where auto-correlations are indeed clearly visible.

Suggestions to implement TE mainly differ in their approaches to remove the
self-predictability of the target signal from the quantification of the
effect. While it has been argued previously that simplistic approaches relying
on a target past operationalisation consisting of a single delay are
insufficient (Wibral et al., 2013), such approaches remain popular. This might
imply that the literature is missing more intuitive demonstrations of the
shortcomings of such approaches. Furthermore, it also has been shown
previously that prominent proposals relying on multidimensional embeddings of
the target time series fail in scenarios of narrow-band effects (Wollstadt et
al., 2017). While suggestions exist that try to overcome this problem by means
of constructing frequency-specific surrogate data (Pinzuti et al., 2020) or
state-space models (Faes et al., 2017), more intuitive explanations of these
failure cases are arguably still lacking. Such an intuitive understanding
should pave the way towards both more widespread awareness as well as simple
fixes of these issues.

One interesting use case for measures that overcome the problems outlined
above is the heavily studied phenomenon of speech envelope tracking as
observed in magneto- and electroencephalography (henceforth MEEG to denote
both modalities) recordings (Ahissar et al., 2001; Hertrich et al., 2012;
Gross et al., 2013; Ding & Simon, 2012; O’Sullivan et al., 2015; Di Liberto et
al., 2015; Wöstmann et al., 2017; Brodbeck et al., 2018; Daube et al., 2019;
Obleser & Kayser, 2019; Zan et al., 2020; Donhauser & Baillet, 2020). In
short, the low-frequency portion of MEEG signals is reliably related to the
time varying energy of the speech signal at a certain delay. In noninvasive
recordings, speech envelope tracking is heavily studied in the canonical delta
and theta bands (Ding & Simon, 2014). At higher frequencies, the effect
usually fails to robustly exceed noise thresholds (but see Kulasingham et al.,
2020). Limiting the analysis to the delta and theta frequency bands is thus an
efficient means to achieve stronger effects. Theories exist that link this
phenomenon to frequency specific mass-neuronal processes of debated
algorithmic significance (Giraud & Poeppel, 2012; Hyafil et al., 2015).

From the perspective of directed connectivity, an interesting question is to
what degree the resemblance of MEEG signals and the speech envelope at a given
time $t$ can be explained from the MEEG signal’s own past, and to what degree
it can only be explained when additionally considering the stimulus. Given the
auto-correlated nature of both neuronal processes and the speech envelope in
the relevant spectral range, it is in theory possible for a system to predict
the upcoming speech envelope (and, as made evident by the progress of past and
current machine learning approaches, also richer parts of the speech signal
Chung et al., 2020; Lakhotia et al., 2021). According to popular theories of
brain function, such predictive coding is also of high utility for a
biological system (Rao & Ballard, 1999; Friston, 2005). The extent to which
speech tracking as quantified by undirected measures such as delayed MI would
decrease when using suitable directed measures would highlight the extent to
which the heavily studied tracking might in fact reflect predictive rather
than reactive processing. This would add to accounts that characterise low-
frequency oscillations as the deliberate effort of biological systems to be in
a state of optimal neuronal excitability (Lakatos et al., 2008; Henry &
Obleser, 2012), such that metabolically costly states of high encoding
fidelity co-occur with the relevant parts of the stimulus (Jones & Boltz,
1989; Schroeder & Lakatos, 2009; Kayser et al., 2015; Młynarski & Hermundstad,
2018).

Here, we consider a simple simulated system of band-limited delayed bivariate
interactions with a clear and intuitively comprehensible ground truth spectral
range and delay. We implement delayed MI, two classic algorithms to estimate
TE as well as a novel, whitening based estimator (“Directed Information based
on conditional entropy”, “DIce”) within the Gaussian copula MI framework
(gcmi, Ince et al., 2017) and extensively test these measures with
simulations. Finally, we explore the behaviour of DIce in an MEG dataset of
continuous speech listening. We find that estimates of the delay between the
stimulus and the response as well as the recovered interaction strength in the
delta and theta bands differ from those recovered by the bivariate delayed MI.

## Results

### An intuitive overview over multiple measures

A first goal of this study was to provide an overview of various
implementations of TE in a simple and intuitively understandable example case.
To do so, we simulated coupled systems with a known ground truth spectral
interaction range and delay. These consisted of 4 - 8 Hz narrow-band filtered
Gaussian white noise to obtain a band-limited source signal, which was delayed
in time by .15 s to obtain a target signal. We added Gaussian white noise to
these signals to mimic the noisy measurement process in neuroimaging. For such
a simplistic but intuitive system, a suitable measure of TE should peak at the
ground truth simulated delay and spectral range. Further, for such highly
auto-correlated signals, TE should yield a highly reduced effect size in
comparison to undirected measures such as delayed MI: a potential source
signal can only add little information if the target signal is already highly
predictable from its own recent history.

Figure 1 shows the results of applying delayed MI, classical TE estimators
TE1D and TESPO as well as our proposal, DIce (see below for more detailed
explanations of each individual measure) to the same time series, both without
and with applying an analysis filter.

A first observation is that in the present implementation with gcmi (Ince et
al., 2017), all measures yield increased effect sizes when a suitable analysis
filter is applied. This, potentially trivially (but see Pinzuti et al., 2020),
demonstrates that the effect size and consequently the sensitivity can benefit
from the application of analysis filters. It thus underscores the utility of
directed connectivity measures that behave robustly when applied to narrow-
band signals. Further, the pass-band of the effect can be found by applying
analysis filters of varying centre frequencies: all measures return the
highest effect sizes across analysis bands within the pass-band of the
simulated effect.

Figure 1: Comparison of delay profiles of various undirected and directed
dependency measures in a simplistic simulation scenario. The simulation here
consists of a 4 - 8 Hz bandpass filtered white noise source signal that is
temporally delayed to obtain the target. White noise is added to both source
and target to model measurement noise. First and middle columns show time
domain sections of the analysed time series, where a filled circle denotes a
target present sample, an empty circle denotes a source past sample at a delay
corresponding to the simulated effect, and crosses denote the respective
target past samples. 2nd and 4th columns show delay profiles of the
respectives measures relative to the ground truth interaction delay and noise
thresholds. 5th column shows spectrotemporal maps, where the frequency
response of the filter used to generate the ground truth effect is overlaid as
a black curve. A Delayed Mutual Information correctly recovers the ground
truth interaction delay. B TE1D also correctly recovers the interaction delay,
but measures an interaction effect that is highly similar to delayed MI. C
TESPO recovers an interaction delay different from the ground truth, even more
so when the analysis is performed on the filtered signal. D DIce only finds a
super-threshold effect with correct recovery of the interaction delay when an
analysis filter is used, but with a strongly reduced effect size in comparison
to other measures.

Secondly, we observe that TE1D fails to return a reduced effect size in
comparison to the undirected delayed MI, even when applied to highly self-
predictable signals as simulated here. This happens for both broadband and
filtered analysis scenarios. TE1D is a classical implementation of TE
(Besserve et al., 2010; Lobier et al., 2014; Ince et al., 2015; Park et al.,
2015; Giordano et al., 2016; Morillon & Baillet, 2017) that operationalises
the target past by means of using only one delay of the target time series
(more specifically, the same delay as that of the source signal relative to
the target present when scanning across delays). For the ground truth delay
simulated here, this single delay fails to capture most of the self-
predictability of the target signal, resulting in an overestimation of the
directed effect.

A third measure, TESPO (Wibral et al., 2013), yields the anticipated strong
reduction in effect size by roughly an order of magnitude in comparison to MI.
It achieves this by using a more effective handle on the target past: a multi-
dimensional “non-uniform” (Vlachos & Kugiumtzis, 2010; Faes et al., 2011)
embedding optimised for self-prediction. As opposed to TE1D, this essentially
consists of multiple delays that are independent of the analysis delay in the
scanning procedure. However, it fails to recover the ground truth interaction
delay simulated in the coupled systems (Wollstadt et al., 2017), in both
broadband and filtered analysis settings. Given this failure in a simplistic
problem setting that is ubiquitous in neuroscientific datasets, it is unclear
to which degree TESPO can practically live up to its promise of correctly
quantifying “predictive information transfer” as well as the interaction
delay.

Finally, the last measure considered in and proposed by this study, DIce,
solves both of these problems: It returns the smallest effect sizes in this
comparison, even to the degree that it fails to detect the effect in a
broadband analysis setting. However, when a suitable analysis filter is
applied, it returns a delay profile with a super-threshold peak at the
simulated delay. It achieves this in a two-step approach (Haugh, 1976): First,
both source and target signals are transformed into time series of surprisal,
i.e. the sample-wise entropy of each time point conditional on the same non-
uniform embedding used in TESPO. Secondly, delayed MI is computed for these
whitened time series.

Taken together, this first analysis provides an overview of four different
connectivity measures in an intuitive and simplistic simulation setting,
illustrating how DIce succeeds in returning intuitive results where TE1D and
TESPO fail.

### Synergy of source and target past about target present distorts
conditional MI based TE

Figure 2: Partial information decomposition perspective on TE. A Delay profile
of DIce for reference. B Delay profile of delayed MI. C Delay profile of
TESPO, which is based on conditional MI and can, according to PID, be seen as
the sum of synergy (F) and unique information of the source (H). D Delay
profile of co-information (abbreviated as “co-info”), which quantifies the
triple set intersection in classic Venn-diagram conceptualisations of three
variable systems. According to PID, it is the net sum of synergy (F) and
redundancy (E). Axis limits chosen to highlight the negative (net synergistic)
portion at positive delays. E - H Delay profiles of the PID atoms redundancy,
synergy and unique information (of target and source). I Delay profiles of
synergy and sign-flipped unique information of the source, highlighting a
surplus of synergy at positive delays. J Delay profiles of DIce when
simulating data with different ground truth effect bandwidths. DIce recovers
the correct delay at all bandwidths (left) and recovers an increasing effect
size as a function of bandwidth. It exceeds the noise threshold in all cases.
K Delay profiles of unique information of the source (left) and synergy
(right) for the same effect bandwidths as in J. For narrow bandwidths, unique
information of the source fails to become positive and recovers the wrong
delay. Further, the sum of unique information of the source and synergy is
dominated by synergy, which features a greater delay estimation error than
unique information.

We wanted to investigate possible explanations for the counterintuitive
results returned by our implementation of TESPO. To do so, we considered a
perspective on TE offered by partial information decomposition (PID). PID is
an information theoretic approach to study trivariate relationships (Williams
& Beer, 2010; Ince, 2017). The central goal in PID is to quantify shared
information (redundancy) between two source variables about a third target
variable. Further, PID aims to measure unique information of both source
variables as well as synergistic information that is only obtainable when
jointly considering both sources. A key insight of PID relevant to TE is that
the basis of TE, conditional MI, is the sum of two “atoms” of PID: unique
information and synergy (James et al., 2016). In other words, conditioning a
bivariate relationship on a third variable will remove redundant information,
but will not deliver solely unique information. It has been pointed out that
this conflation of unique and synergistic information defies interpretations
of TE as measuring “information flow” (James et al., 2016). By relying on
conditional MI, TE measures not just unique information of the source about
the target present ignoring the target past, but instead conflates this with
synergistic information stemming from interactions of the source and the
target past. The resulting quantity is thus not “localisable” to the source
(James et al., 2016).

We turned to data from the same simulated system as analysed in the previous
section. In a first step, we considered co-information (McGill, 1954) of the
source, the target present and the target past. It can be obtained by
subtracting conditional MI from bivariate MI, and thus quantifies the triple
set intersection in classic Venn diagram conceptualisations. It can take on
positive and negative values, where positive values denote redundancy and
negative values denote synergy. It is important to note that from a PID
perspective, the conceptually simpler (and less controversial) co-information
conflates redundancy and synergy to a single net quantity, which PID aims to
decompose into pure redundancy and synergy. Figure 2D shows negative co-
information values at later delays, demonstrating that there is indeed a net
synergy of the source and the target past about the target present.

Next, we applied PID using an implementation based on common change in
surprisal (Iccs Ince, 2017). We find that TESPO (figure 2C), as it is based on
conditional MI, can be decomposed into largely identical profiles of synergy
(figure 2F) and unique information of the source (figure 2H) differing only in
their sign. The positive net conditional MI (i.e. TESPO) however turns out to
entirely stem from a surplus of synergy at later delays that is not cancelled
out by unique information of the source (figure 2I). In other words, in this
simplistic example case, TESPO measures a synergistic effect.

As a consequence of the PID perspective on conditional MI based TE
implementations, it has been proposed to use the unique information of the
source as a more appropriate measure in order to avoid the quantification of
synergistic effects (Barrett, 2015). However, in our example case, this
quantity is negative across the entire range of considered interaction delays
(figure 2H), meaning that from considering the source on its own, predictions
of the target would become worse. This can be seen as the result of two
factors: firstly, in this simulation, the source signal is noisy, and
secondly, the simulated effect has a very narrow spectral range and thus
highly limited degrees of freedom. As these factors come together, the
efficient operationalisation of the target past makes it impossible to improve
on its prediction of the target present.

We thus considered similar simulations with varying bandwidths of the ground
truth effect, and verified that at broader bandwidths, the unique information
of the source becomes positive (figure 2K). Of note, the delay recovered by
unique information of the source generally has lower deviations from the
ground truth than what is recovered by synergy, but this deviation is still
non-zero at narrower effect bandwidths. We further find that with increasing
bandwidths of the effect, synergy increases as well and peaks at delays closer
to the ground truth delay. Crucially however, across all simulated bandwidths,
our proposed measure, DIce, finds super-threshold effects at the correct delay
(figure 2J).

Taken together, we have shown that the delay mis-estimation problems of TESPO,
when applied to narrow-band signals, stem from the known conflation of unique
and synergistic effects in conditional MI based TE implementations. In this
example, the effect quantified by TESPO is in fact dominated by synergy. We
further show that a possible alternative, namely unique information of the
source as obtained from PID, has both a lower sensitivity and a worse
performance in recovering the ground truth interaction delay than our proposed
conditional MI-free measure DIce.

### Varying the ground truth interaction delay

After these initial intuitive and single sample based demonstrations of
problems of two classic algorithms to quantify TE as well as a suggestion to
explain their origin, we were interested in more thoroughly testing the
characteristics of the measures. To do so, we performed extensive analyses of
the simplistic coupled systems that measured the problems over repeated
samples. For each of these samples, we computed the delay profiles of the set
of four measures and evaluated the recovered interaction delay as well as the
recovered effect size as obtained from the peaks across the delay profile of a
given repetition.

Figure 3: Simulation with varying ground truth delay. Basic setup of the
simulation is the same as in figure 1, however, here, the ground truth
interaction delay is varied, and 100 repetitions are sampled. Moreover, a
higher SNR is used. Plots show the median across repetitions, shaded regions
indicate bootstrapped 95% confidence intervals. A Recovered effect size when
the data are analysed without a filter. TE1D suffers from recovering a
systematically varying effect size across different ground truth delays
despite a constant simulated interaction strength. B Same as A, but applying
an analysis filter. C Recovered delays when data are analysed without a
filter. TESPO underestimates the true interaction delay. D Same as C, but
applying an analysis filter.

In a first simulation, we were interested in the characteristics of the
measures across a range of simulated ground truth interaction delays. In
principle, a suitable measure of TE should recover the same constant effect
size at varying simulated interaction delays when all other factors are kept
constant. In figure 3, the results demonstrate that this is indeed the case
for all measures in both broadband and filtered analysis scenarios except for
TE1D. This measure exhibits a systematic variation in the recovered effect
size across different simulated interaction delays reminiscent of the filter
ringing of the auto-correlation profile of the target signal. The failure here
again highlights the problems of operationalising the target past with a
single delay that varies as the analysis delay is scanned (Wibral et al.,
2013). In case of very short ground truth delays, the MI at the peak of the
delay profile is conditioned on the same short delay of the target variable.
Since narrow-band signals have high auto-correlation at such short delays,
this leads to a strong reduction of conditional MI vs MI. As the ground truth
delay increases, the MI is conditioned on longer delays, where the auto-
correlation of the target variable wanes and waxes and thus leaves a
correspondingly varying conditional MI. This can lead to potential
interpretational pitfalls when for example the results of TE1D obtained in two
different experimental conditions are compared. If these conditions simply
differ in delay, this will lead to different recovered effect sizes and could
thus be falsely interpreted as differences in directed dependency.

Further, the results reiterate that TESPO, while featuring a constant effect
across different ground truth interaction delays, recovers a flawed estimate
of the interaction delay. DIce on the other hand behaves favourably by
returning a constant effect size at the correct interaction delay. However,
due to its drastically reduced effect size in this band-limited interaction,
its estimates are noisier.

### Varying the signal-to-noise ratio in source and target signals

Figure 4: Simulation with varying signal-to-noise ratios. Basic setup of the
simulation is the same as in figure 1, however, here, the amplitude of noise
added to source and target varies (but the same amount is added to source and
target), and 100 repetitions are sampled. Plots show the median across 100
repetitions, shaded regions denote bootstrapped 95% confidence intervals. A In
a broadband analysis scenario, all measures recover decreasing effect sizes
with increasing noise amplitudes. DIce recovers the smallest effect sizes and
reaches the noise floor the earliest. B In a filtered analysis scenario, the
recovered effect sizes are in general higher for all measures. C The recovered
delays deteriorate as the noise increases. In a broadband analysis scenario,
all measures correctly recover the ground truth interaction delay at low noise
amplitudes. D In a filtered analysis scenario, TESPO has higher delay
reconstruction errors, while the delays recovered by the other measures
benefit from the increased effect size.

In a second simulation, we were interested in comparing the four measures when
faced with increasingly noisy signals. While the ideal measure should be
highly sensitive to a present effect and recover it at the correct delay,
ignoring self-predictable parts of dependencies should inevitably reduce the
effect size.

We indeed found that DIce recovered the smallest effects and hit the noise
floor at the lowest noise amplitude in comparison to the other measures
(figure 4). In a filtered analysis scenario, the recovered effects were in
general higher than those in a broadband analysis scenario. Interestingly, for
some noise amplitudes, DIce recovered stronger effects than TESPO. At very low
noise amplitudes, all measures succeeded in recovering the correct delay
irrespective of whether filters were applied in the analysis or not. However,
as noise increased, especially TESPO failed in recovering correct delay
estimates. While this was especially problematic for filtered analyses, other
measures had lower delay estimation errors at higher noise amplitudes due to
the gain in sensitivity afforded by filtering.

### Varying the signal-to-noise ratio independently in source and target
signals

In a third step, we were interested in assessing how the four measures would
react to asymmetric variations of the SNR (Bastos & Schoffelen, 2016). The
ideal measure of TE should recover an invariant interaction delay and an
effect size that symmetrically decreases as noise of increasing amplitude is
added to either source or target signals.

Figure 5: Simulation with signal-to-noise ratios varying independently in
source and target signals. Basic setup of the simulation is the same as in
figure 1, however, here, the amount of noise added to source or target signals
is varied, and 100 repetitions are sampled. A - D show single trial delay
profiles in a broadband analysis scenario for various measures, corrected for
the true interaction delay across different noise amplitudes (negative delta
on the y-axis refers to higher noise amplitude in the target signal, positive
delta on the y-axis refers to a higher noise amplitude in the source signal).
The ideal measure should always return the highest values at a lag of 0
relative to the true delay. Only MI and TE1D succeed (but do not quantify
Granger-causal information). E and F show the same as A – D, but repeated for
100 trials. Plots show the median across trials, shaded regions denote
bootstrapped 95% confidence intervals. The maximum effect size across lags for
all measures symmetrically decays as the noise increases in either source or
target signals. The delay misestimation however behaves asymmetrically,
especially for TESPO, which recovers wrong delays when the source is noisier
than the target. G – J The same as A – D, but within a filtered analysis
scenario. Here, DIce performs favourably. K and L show the same as E and F,
but within a filtered analysis scenario.

Figure 5 shows how the set of four measures behaved. All measures succeeded in
returning symmetric decreases of the recovered effect sizes when there was
more noise in either source or target signals, irrespective of whether
analysis filters were applied or not. However, only the undirected delayed MI
as well as the practically undirected TE1D (see figure 1) recovered
interaction delays that were unaffected in these settings. For TESPO, an
increasing SNR imbalance with noisier source signals led to stronger biases in
the estimated interaction delay for broadband analyses. For analyses where
filters were applied, these biases grew stronger for both noisier source and
noisier target signals. DIce exhibited biases in the broadband analyses, but
performed favourably when analyses filters were applied. This analysis
additionally corroborates an argument of caution when interpreting TE results
(Wollstadt et al., 2017): when there are asymmetric changes in e.g.
measurement noise across conditions, this does not automatically imply a
change of coupling across these conditions, but might merely reflect the
sensitivity of TE to such measurement noise.

### Varying the bandwidths of the simulated effect as well as of the analysis
filter

In our fourth simulation, we reasoned that a characteristic effect which a
suitable measure of TE should show is a varying recovered effect size as a
function of the bandwidth of the ground truth effect. Specifically, the
theoretical extreme case of a sinusoidal source signal and its phase-shifted
copy as a target signal should lead to zero TE. Such signals have no degrees
of freedom and are perfectly predictable from themselves. Additionally
considering potential source signals can thus not add any information. As one
turns to signals with increasingly broad passbands, these degrees of freedom
increase. Consequently, a higher directed effect should be quantifiable.

Figure 6: Simulation with varying bandwidths of the transmitted signal. Basic
setup of the simulation is the same as in figure 1, however, here, the ground
truth bandwidth of the effect is varied, and 100 repetitions are sampled.
Plots show the median across 100 repetitions, shaded regions denote
bootstrapped 95% confidence intervals. A The recovered effect size increases
as the bandwidth of the effect grows. This increase however spans the most
orders of magnitude for DIce. B In a filtered analysis scenario, the effect
sizes are stronger than in a broadband analysis scenario. C When varying the
bandwidth of the filter used in the analysis, DIce returns the peak effect
size when the analysis filter matches the ground truth effect. All other
measures fail to do so, TESPO even increases monotonically as the analysis
filter passband grows. D - F report the interaction delays recovered in the
same simulations as shown in A - C. The previously reported misestimation of
the interaction delay of TESPO is limited to narrow band effects.

When testing this with the set of four measures (figure 6), we found that all
measures monotonically increased in effect size as the bandwidth of the effect
increased. However, in a broadband analysis, this increase spanned 2.5 orders
of magnitude for DIce, which thus had the strongest relationship between its
recovered effect size and increased degrees of freedom of the input signals.
This was however also driven by the ground truth effect covering a larger part
of the entire frequency spectrum. In a filtered analysis scenario where filter
parameters were chosen to isolate the effect from the noise, only TESPO and
DIce showed increases in effect size with broader effect bandwidths. The
recovered delays were, as shown in previous analyses, unaffected for all
measures except for TESPO. The latter exhibited stronger biases in the
recovered delay for narrower simulated effect bandwidths.

In theory, applying analysis filters that match the ground truth effect should
return a stronger effect size than analysis filters that are too broad or too
narrow. Too broad analysis filters should fail to reduce the impact of noise
in neighbouring frequencies, while too narrow analysis filters should ignore
variance of the effect of interest and thus reduce the sensitivity. When
applying analysis filters of varying width to a simulated coupled system of
fixed ground truth effect bandwidth, DIce showed the clearest peak close to
the ground truth effect (bandwidth of 10 Hz, figure 6C) across analysis filter
widths. TESPO on the other hand exhibited a monotonic increase in effect size
as the analysis filter was broadened, which would thus result in mis-
estimations of the interaction bandwidth. The undirected delayed MI and the
practically undirected TE1D only showed comparably weak variance across
different analysis filters.

### Varying the centre frequency of the simulated effect

In a last simulation, we were interested in assessing the stability of the set
of four measures when the centre frequency was varied. A suitable measure of
TE should recover effect sizes and interaction delays with no variation under
different ground truth centre frequencies, and it should find the strongest
effects at the ground truth frequency.

Figure 7: Simulation with varying centre frequencies of the effect. Basic
setup of the simulation is the same as in figure 1, however, here, the centre
frequency of the effect is varied, and 100 repetitions are sampled. Line plots
show median across 100 repetitions, shaded regions denote bootstrapped 95%
confidence intervals. A In a broadband analysis scenario, all measures recover
effect sizes that are highly constant across centre frequencies. TE1D however
exhibits a slight ringing, TESPO and DIce recover noisier estimates. B The
same holds for a filtered analysis scenario, where DIce and TESPO exhibit
stronger variations. C and D As shown in other simulations, TESPO suffers from
delay estimation problems. These vary especially strong across centre
frequencies for filtered analysis scenarios. E - H Confusion matrices obtained
when scanning data with different ground truth interaction frequencies with a
bank of analysis filters (choosing the frequency with the maximum across
analysis filters). All measures generally succeed in returning peak effects at
the true interaction frequency. TESPO has the highest error.

We found that all measures considered here recovered effect sizes that were
indeed relatively stable across different ground truth centre frequencies in
both broadband and filtered analysis scenarios (figure 7). TE1D and DIce
however did exhibit slight biases, especially for the combination of DIce and
higher simulated centre frequencies. All measures except TESPO recovered
unbiased estimates of the interaction delay for all tested centre frequencies.
TESPO on the other hand exhibited a marked variation of the recovered
interaction delay across the different centre frequencies.

Finally, when analysing the simulated data with filters across different
bands, all measures mostly peaked when analysis filters matched the ground
truth frequency. For MI and TE1D, the mean absolute error (MAE) across
frequency bands was 0 frequency bins, for DIce it was .01 frequency bins and
for TESPO, the MAE was strongest with .04 frequency bins. Identifying the
correct spectral band thus did not seem to pose a particular problem for any
of the considered measures.

### Studying low-frequency MEG speech envelope tracking with directed measures

Finally, we wanted to explore the behaviour of our measure on real data. To do
so, we turned to a dataset of $n=24$ participants who listened to a continuous
narrative of $1$ hour duration while their MEG was recorded (Daube et al.,
2019). Examples of cerebro-peripheral coupling such as this present a good
testbed for frequency-specific measures of directed connectivity, because the
ground truth direction of the effect is known. Further, a plethora of studies
has examined the relationship of MEEG responses to the time varying energy, or
“amplitude envelope”, of the speech stimulus in the delta and theta bands
(Ding & Simon, 2014). However, it is usually studied using bivariate,
undirected measures of connectivity that do not consider the self-
predictability of the MEG responses. The degree to which such bivariate
dependencies could be accounted for by auto-regressive models of the neuronal
response signal could in principle reflect the degree to which bivariate
speech tracking is a signature of predictive rather than reactive bottom-up
processing. On the other hand, effects found by directed measures of
dependency (under the assumption of a given auto-regressive model) would be
stronger evidence of reactive, bottom-up processing of unpredictable parts of
the stimulus input.

Figure 8: Caption on following page. Figure 9: Results obtained from source
level MEG recordings (left and right auditory cortices, l and r ACs) obtained
during continuous speech listening (n=24). A Results for a typical
participant. Plots show spectrotemporally resolved delayed MI with individual
boundaries of delta and theta bands overlaid as black horizontal lines
(leftmost plot) as well as delay profile of delayed MI and DIce in delta and
theta bands, with noise thresholds overlaid as black dotted lines. B Spectral
profiles of delayed MI for all participants in left (top) and right (bottom)
ACs. Points denote upper pass-band cut-off frequencies of the delta (upper
scatter) and theta (lower scatter) bands. C Binarised delayed MI as well as
DIce delay profiles for all participants in delta and theta bands in left and
right ACs. While delayed MI in both bands and DIce in the theta bands cross
the noise threshold in most cases, DIce in the delta band is often close to or
below the noise threshold. D Ratios of DI/MI in delta and theta bands of left
and right ACs for each individual participant, medians overlaid as black
lines, density estimates from posterior distributions overlaid as transparent
shapes. Theta band bivariate tracking can be less explained by auto-prediction
of the MEG signal than delta band tracking. E Effect sizes underlying the
computation of the ratio in D in delta and theta bands of the left and right
ACs in each individual participant. Medians are overlaid as coloured lines,
density estimates from posterior distributions overlaid as transparent shapes.
F Delays recovered by delayed MI and DIce in the delta and theta bands of left
and right ACs. Since DIce fails to cross the noise threshold in many cases for
the delta band, these delay estimates are uninterpretable. In the theta band,
DIce recovers a longer delay than delayed MI. Also see Figures S1 – S4.

In most participants, we found spectra of delayed MI (see figure 9B and figure
S1) that were suggestive of two spectral components involved in speech
tracking which we will refer to as delta and theta bands (note that our
functional definition of the theta band had higher upper cutoffs than the
canonical theta band, which is commonly defined between 4 to 8 Hz, Klimesch,
1999; Wang, 2010). We found super-threshold delayed MI in both left and right
auditory cortices of all participants in both delta and theta bands (figure
9C). This was stronger in delta than in theta frequency bands (fraction of
samples in favour of hypothesis $f_{h_{1}}=1$, figure 9E, see also figure S4).
With DIce however, we found only weak effects that barely exceeded the noise
thresholds in the delta bands, while effects in the theta bands exceeded noise
thresholds with only 1 exception in the left and 2 exceptions in the right
hemisphere (figure 9C). This constitutes strong evidence for a robust
population-level prevalence of the effect (Ince et al., 2021). These theta
DIce effects were generally much stronger than in the delta band
($f_{h_{1}}=1$; figure 9E, see also figure 9E, see also figure S4). These
results translated into a difference of ratios of DI divided by MI between the
delta and theta frequency bands, where we found higher ratios in the theta
than in the delta band in both left and right auditory cortices
($f_{h_{1}}=1$; figure 9D, see also figure 9E, see also figure S4). This
suggests that the bivariate speech tracking in the theta band could consist of
more bottom-up and reactive processing than the bivariate speech tracking in
the delta band. We had found in simulation analyses that the effect size
recovered by DIce grows as a function of increasing bandwidth (see figure 6).
Since our individualised definitions of delta and theta bands led to wider
theta than delta analysis filters, we repeated this analysis while
constraining the analysis filter bandwidth to be the same in delta and theta
bands. We found that the difference in DI over MI ratios of the delta vs the
theta bands shrunk, but persisted under this constraint ($f_{h_{1}}=1$; see
figures S3 and S4). We could thus rule out that analysis bandwidth alone could
explain the difference in DI/MI ratios in delta vs theta bands.

Given the overall relatively narrow passbands of both delta and theta bands,
the measurable directed effects were in general up to several orders of
magnitude lower than the undirected effects, suggesting that predictive
processing could indeed make up the lion’s share of delta and theta envelope
tracking.

Lastly, we considered the delays recovered by delayed MI and DIce in delta and
theta bands (figure 9F). For the delta band, the delays recovered by DIce
spanned a wide range, even reaching into regions suggestive of the MEG signal
preceding the speech signal for some participants. However, since the effect
sizes in the delta band were in many cases close to or below the noise
threshold, these recovered delays were uninterpretable. For the theta band
however, effect sizes were robust in most cases. Interestingly, the delays
recovered by DIce (lAC: median across participants of -.43s, rAC: median of
-.40s) were longer than those recovered by delayed MI (lAC: median of -.10s,
rAC: median of -.10s; $f_{h_{1}}=1$), suggesting a slower bottom-up, reactive
processing.

Taken together, we take our results as aligning with models that characterise
the lion’s share of activity of auditory cortices visible in MEG at a given
point in time as reflecting predictions of the stimulus at a short delay. Such
predictions can be formed on the basis of an integration of reactive, bottom-
up processing of stimulus input at longer delays with an internal language
model.

## Discussion

In this study, we have addressed the problem of quantifying band-limited
directed interactions in bivariate sets of source and target signals. With a
simplistic and intuitive simulation setup, we have highlighted shortcomings of
common estimators of TE when facing this problem. Our proposed alternative
DIce, relying on a type of temporal whitening of the source and target
signals, overcomes these shortcomings. With an array of simulations, we
extensively characterised DIce in relation to undirected MI as well as common
estimators of TE. Lastly, we turned to a dataset of continuous speech
listening in MEG to study the speech envelope tracking specific to delta and
theta bands from the perspective of a directed measure. We found that in such
narrow bands with essentially low degrees of freedom, the measurable directed
effects are very small in comparison to undirected effects. Moreover, with
DIce, we found that theta speech tracking has a stronger directed part than
delta speech tracking, and that the directed theta speech tracking has a
longer delay and a broader temporal profile than the undirected tracking. The
directed effects could potentially reflect a purely bottom-up processing stage
at a later delay. On the basis of this reactive processing, the auditory
cortical system could then generate predictions about the upcoming speech
acoustics, which the undirected effects could be a signature of.

Many of the problems and perspectives highlighted in this study are not new to
the literature. We are for example not the first to point out problems of TE
estimators relying on only a single sample of the target past signal at a
delay equivalent to that of the source to the target at a given scanned
interaction delay (Wibral et al., 2013). However, we hope that our intuitive
demonstrations contribute to a more widespread appreciation of the high
similarity of the resulting “TE”1D and simple delayed MI as well as variations
of recovered effect sizes in the presence of different ground truth
interaction delays. Likewise, the delay misestimation problems of more
sophisticated estimators such as TESPO when faced with narrow-band effects
have been reported before (Wollstadt et al., 2017). In the same way, the PID
perspective on conditional MI based TE and the implication of including
synergistic interactions has been developed previously (James et al., 2016).
Here, we found that from this PID perspective on the TE problem, the delay
mis-estimations of TESPO could mostly be attributed to such synergistic
interactions of the source and the target that TE estimators relying on
conditional MI pick up on.

An intuitive solution to fix this would be to instead consider the unique
information of the source about the target, and thus effectively ignore the
synergistic contributions included in conditional MI (Barrett, 2015). However,
we found that this was neither particularly sensitive in detecting effects nor
accurate in recovering the ground truth delays, albeit with much less error
than synergy-driven conditional MI based estimators. Our proposed estimator,
DIce, circumnavigates these issues by not relying on conditional MI in the
first place, but by instead converting the source and target time series into
time series of sample wise conditional entropy. This effectively frees the
signals from predictable parts, ensuring that a subsequently computed delayed
MI cannot be affected by autocorrelation. Applying this whitening to the
target signal only would be analogous to the basic idea of TE, but would
introduce asymmetric temporal distortions in the target relative to the
source, resulting in erroneous delay estimates. Again, the general idea of
such a two-staged approach is not new as such (Haugh, 1976; Cliff et al.,
2021). However, we have here developed it from a simple cross-correlation of
residuals (Haugh, 1976), which only takes into account the mean of a predicted
time series, to the delayed MI of sample wise conditional entropy, which
instead relates a given sample to the mean and variance of a predicted
distribution. On the basis of our simulations, we put forward that it does
offer a useful perspective on the problem of estimating directed interactions
of band-limited processes at the correct delay.

In the light of various careful considerations of the pitfalls associated with
this goal (Florin et al., 2010; Barnett & Seth, 2011), we thus hope that we
can contribute to a more widespread adoption of directed perspectives on band-
limited dependencies. Within neuroscience, such genuine narrow-band signal
components are arguably not in all places where aperiodic signals have been
straight-jacketed into oscillations by means of band-pass filters (Donoghue et
al., 2020; Gerster et al., 2021). In fact, depending on the algorithm used to
define oscillatory components, even the delta component identified by us is
potentially attributable to the aperiodic part of the MEG spectrum. In
general, oscillatory components are nevertheless indisputably widespread
(Wang, 2010). It is in such situations of high auto-correlation where the
consideration of self-predictability should have the most obvious and
strongest effects, and where the use of measures like TE is thus of
potentially high interest. After all, frequency-specific oscillatory
components are at the heart of popular conceptualisations of interactions of
neuronal populations (Schnitzler & Gross, 2005; Fries, 2015). We conjecture
that estimators of directed interactions that return counterintuitive results
in such circumstances as simplistically simulated may not be well matched to
the interpretations typically applied. This problem has been recognised
elsewhere, and (Pinzuti et al., 2020) propose a clever permutation scheme to
address this issue. Here, we instead suggest to prioritise the
interpretability of the directly measured effect size.

In this spirit, our measure is in principle a versatile tool applicable to
arbitrary (neuro-)scientific questions. We however subscribe to the view that
as such, mechanistic interpretations of cerebro-cerebral dependencies are
prohibitively hard in many cases (Mehler & Kording, 2018). Given the highly
incomplete picture of the entire neuronal activity that is accessible with any
given neuroimaging modality, it is virtually impossible to rule out that a
given dependency stems from a third unobservable region affecting the supposed
source and target. Directed dependency measures per se thus never warrant the
inference of causality. This problem is however alleviated in cerebro-
peripheral settings, where the dependency of neuroimaging signals on external
signals is studied (Gross et al., 2021).

The question of causality is thus arguably less controversial in our example
case of the dependence of MEG signals on a continuous acoustic speech
stimulus. According to the burgeoning field of predictive coding however, the
brain appears to predict the upcoming speech input (Brodbeck et al., 2018;
Donhauser & Baillet, 2020; Heilbron et al., 2021), rendering temporal
relations of time series of assumed cause and effect less trivial. Our
suggestions regarding these problems are mostly of an indirect nature, arguing
that the part of undirected dependencies that directed measures deliberately
ignore could reflect predictive processing, while what directed measures
quantify should reflect reactive, bottom-up processing. It is important to
note however that simply because the part of a neural response that we can
observe with a given neuroimaging modality is predictable from its own past
does not automatically imply that it is in fact being predicted by the brain
(de-Wit et al., 2016). Response components identified with impulse response
functions as measured with autocorrelations, temporal response functions
(Crosse et al., 2016) or here delayed MI have been suggested to be analogous
to the well-studied components of auditory event-related potentials (Lalor et
al., 2009). These do occur in response to unpredictable events (Ritter et al.,
1968), implying that the underlying generators are not at all restricted to
prediction. It has however also been suggested that the low-frequency
potential “entrains” in response to sustained stimuli such as continuous
speech, especially when they are (quasi-)rhythmic (Lakatos et al., 2008,
2019). Viewing ERPs in response to unpredictable events through this lens,
they could manifest the launch or “reset” of such predictions for subsequent
events (Sayers et al., 1974; Makeig et al., 2002). Such entrainment would
serve to align phases of low and high neuronal excitability to relevant parts
of the stimulus, such that crucial stimulus information could be encoded with
high fidelity. In this sense, the low-frequency potential is indeed often seen
as a neuronal signature of a prediction about the upcoming stimulus. A further
interesting detail of our results in this respect is the distribution of lags
recovered by delayed MI (i.e. the timing of the peak of the delay profile).
This was relatively broad in the delta frequency band, with the delays of some
participants even suggesting the stimulus to follow the response, a hallmark
of prediction that has been observed and interpreted similarly before in EEG
delta-band responses to continuous speech (Etard & Reichenbach, 2019).
Interestingly, such slow delta tracking signatures disappear when speech of a
language foreign to the participant is listened to (Ding et al., 2016), which
is to be expected if it is seen as a predictive component relying on an
internal language model.

Under this perspective, the reactive, bottom-up part of the coupling carved
out by DIce would then essentially be driven by adjustments of the predictions
in reaction to unpredictable parts in the stimulus. An important constraint of
our proposed approach however follows from the choice of the model used to
whiten the signals. Our choice of a non-uniform embedding to predict a
signal’s present state from its own past is rooted in the history of TE, and
its suitability to model an actual hypothesised biological acoustic model is
disputable. In principle, there is nothing that prohibits the combination of
our framework with more powerful autoregressive models such as recurrent
neural networks. That the auditory system is actually entertaining
correspondingly complex multi-level predictions of the upcoming input is an
increasingly popular perspective (Donhauser & Baillet, 2020; Koskinen et al.,
2020; Heilbron et al., 2021; Jain et al., 2021; Schmitt et al., 2021;
Caucheteux et al., 2021). This would then effectively further shrink the
effect sizes that our proposed DIce could discover, assigning an even smaller
portion of variance to unambiguously bottom-up reactive processing than found
here.

Lastly, it is interesting to compare our approach to a recently increasingly
popular analysis strategy in the field of predictive coding. This strategy
usually builds on the framework of encoding models (Naselaris et al., 2011),
where first a so-called linearising feature space (a nonlinear transformation
of the stimulus) is identified and then linearly mapped onto the brain
response. For questions of predictive coding, researchers make use of powerful
predictive models of a stimulus class to derive measures of surprisal (and
uncertainty) associated with a given part of the stimulus given its preceding
context to obtain the linearising feature space (Brodbeck et al., 2018;
Donhauser & Baillet, 2020; Koskinen et al., 2020). Our transformation of the
signal of time varying energy of the speech stimulus into a time-series of
conditional entropy can be seen as a simple version thereof. Within the
standard approach, it is then however uncommon to perform a similar operation
on the neuronal response, which is central to our proposed measure. The
standard approach makes sense from the viewpoint of the hypothesis that the
neuroimaging signal mainly reflects a prediction error. The conceptualisation
of neuronal signals as reflecting prediction errors however decreases in its
appeal to the degree to which the neuronal signal is predictable from itself.
How could the brain be surprised if it already knew it was going to be
surprised? We thus hypothesise that studies of prediction errors could
increase their specificity by including a whitening of the target signals akin
to what we suggest here. In theory, this should allow a more detailed
characterisation of processes related to prediction errors.

Taken together, this study makes a threefold contribution to the problem of
quantifying band-limited directed dependencies: We offer simplistic yet
intuitive and accessible simulations, propose our own estimator based on an
information theoretic pre-whitening as well as provide an application to the
case of speech tracking in MEG in delta and theta bands. With this, we hope to
spark further interest in directed perspectives on frequency-specific
interactions.

## Methods

### Estimation of information theoretic quantities

We used the Gaussian Copula Mutual Information framework (gcmi, Ince et al.,
2017). Here, the basic idea is to transform variables into standard normals to
then apply closed-form expressions for information theoretic quantities of
Gaussians. These can be derived as follows. The entropy $H$ of a variable $X$
is given by the expected value of the surprisal:

$H(X)=E(h(x))$ (1)

For a Gaussian continuous variable $X$, this corresponds to the “global”
differential entropy:

$H(X)=\int_{-\infty}^{\infty}f(x)h(x)dx$ (2)

Here, the local or “sample wise” contribution, i.e. the surprisal (also called
“information content”, “Shannon information” or “self-information”), is thus
given by:

$h(x)=-\log(f(x))$ (3)

For a Gaussian continuous variable of $k$ dimensions, this is equal to to the
negative log-likelihood:

$h(x)=-\left(-\frac{1}{2}\left[\log(|\Sigma|)+(x-\mu)^{\prime}\Sigma^{-1}(x-\mu)+k*\log(2\pi)\right]\right)$
(4)

$H(X)$ can then be simplified to:

$H(X)=\frac{k}{2}+\frac{k}{2}\log(2\pi)+\frac{1}{2}\log(|\Sigma|)$ (5)

We can then compute the mutual information between the variables $X$ and $Y$
(where $[]$ denotes concatenation to obtain a joint variable) follows:

$MI(X,Y)=H(X)+H(Y)-H([XY])$ (6)

This pragmatic estimator comes at the cost of not being able to quantify all
possible nonlinear effects, as information theoretic measures should do in
theory. Instead, it quantifies relationships described by a Gaussian copula.
However, it avoids the loss of information and incompatibility with
multidimensional variables as well as higher-order information theoretic
quantities incurred by binning, a similarly pragmatic choice. Further, it is
computationally cheaper and has a higher sensitivity for rank relationships
than more sophisticated estimators such as nearest neighbour approaches, which
can resolve nonlinear effects, but can in practice only be computed on vast
amounts of samples from variables of relatively limited dimensionality. In
situations where nonlinear effects are of interest and resolvable when using a
binned estimator with a realistically low amount of bins (such that each bin
of the joint distribution is sufficiently sampled), such estimators are a
pragmatic complement to gcmi.

### Delayed mutual information

The simplest information theoretic approach of estimating a dependency between
source and target variables is computing the delayed mutual information (MI).
It is the information theoretic equivalent of cross-correlation, where
versions of the source- and target variables with a range of lags to each
other are generated to then compute the dependency measure at each lag and
obtain a delay profile of the relationship. In principle, the lagging of the
source and target variables to each other can be implemented either way, by
lagging the source and fixing the target or vice versa. To simplify concrete
explanations that follow, we here choose an implementation where the source is
lagged while the target is kept fixed. We thus define delayed MI between time
series $X$ and $Y$ at an interaction delay $\delta_{a}$ as:

$\text{ delayed }MI(X,Y)_{\delta_{a}}=MI\left(X_{\delta_{a}},Y_{0}\right)$ (7)

Here, $X_{\delta_{a}}$ denotes the lagged source with negative values of
$\delta_{a}$ corresponding to delays where the source precedes the target and
positive values corresponding to delays where the source follows the target.
We trim the target present $Y_{0}$ to be of an equal number of samples as the
lagged source.

### Transfer entropy

Originally, transfer entropy (TE) has been described as a conditional MI
(Schreiber, 2000). Instead of computing the simple delayed MI at each lag, the
idea is to condition this MI on an operationalisation of the target past. In
principle, the resulting quantity is thereby supposed to quantify to what
extent the target signal can be explained by the source signal over and above
to the extent to which the target signal can be explained by its own past.
This idea is equivalent to that of Granger causality (Granger, 1969) in the
case of linear effects observed in Gaussian variables (Barnett et al., 2009).
Transfer Entropy based on conditional MI can be expressed as a difference of
two MI terms:

$TE(X,Y)_{\delta_{a}}=MI\left(X_{\delta_{a}},\left[Y_{-}Y_{0}\right]\right)-MI\left(X_{\delta_{a}},Y_{-}\right)$
(8)

A critical problem is the concrete operationalisation of the target past
$Y_{-}$.

#### One-dimensional embedding

A simple form of operationalising the target past is to simply mirror the
delay of the delay scanning procedure of the delayed MI (Besserve et al.,
2010; Lobier et al., 2014; Park et al., 2015; Ince et al., 2015; Giordano et
al., 2016; Morillon & Baillet, 2017). When a given interaction delay is
considered in diagonal TE, the delayed MI between $X$ and $Y$ is conditioned
on a version of the target to which the same lag is applied as that which is
applied to the source. We refer to this one-dimensional embedding as “TE1D”.
With $Y_{-}$ thus corresponding to $Y_{\delta_{a}}$, we obtain:

$TE_{\operatorname{1D}}(X,Y)_{\delta_{a}}=MI\left(X_{\delta_{a}},\left[Y_{\delta_{a}}Y_{0}\right]\right)-MI\left(X_{\delta_{a}},Y_{\delta_{a}}\right)$
(9)

#### Multidimensional embedding approaches

The one-dimensional method of operationalising the target past is simple to
implement, but hard to motivate. When considering that the core goal of
quantifying TE is to estimate to what degree the source variable contains
information about the target that is not available from the target past
itself, it becomes obvious that the assumption of capturing the target’s
autoinformation with a single delay that varies as a function of the
considered interaction delay is daunting (Wibral et al., 2013). Multiple
target past time points could jointly influence the target present, and there
is no obvious reason why the target past should be varied as a function of the
interaction delay considered. A source variable could have information about
the target variable at a different delay than the delay at which the target
past informs its own future in the same way as the source. A more sensible
approach to describe the target past is thus to consider a multidimensional
representation or “embedding” (Takens, 1981) that is invariant to the
interaction delay considered between source and target in the delay scanning
procedure. A classic idea to find such an embedding of the target past is
parameterised by two parameters: A spacing parameter and a dimensionality
parameter (Ragwitz & Kantz, 2002). This means that the assumption is that the
information the target past contains about its present can be captured with a
number of equally spaced delays relative to the present. A further development
of embedding methods relaxes the assumption that the embedding delays must be
uniformly spaced (Vlachos & Kugiumtzis, 2010; Faes et al., 2011). It is
conceivable that a lot of information about the target present can be found in
a combination of many, densely clustered delays at very short time scales as
well as only few unevenly spaced delays at longer time scales, or vice versa.
To account for this, non-uniform embedding procedures adopt an iterative
approach to generate a target past embedding. Concretely, a search space of
candidate delays is defined by the maximal delay that is assumed to have an
influence on the target present. In the first iteration, the delay within this
search space is chosen that maximises the MI of the target variable at this
delay relative to the target present. In subsequent iterations, a conditional
MI between all remaining candidate delays and the target present conditioned
on the already chosen delays is computed. If desired, it can be tested whether
the found maximal conditional MI surpasses a defined noise threshold. In that
case, the iterative procedure can be stopped once the found maximal
conditional MI no longer exceeds the threshold. Alternatively, a pragmatic
hyperparameter of maximal iterations can be defined, which avoids the
computationally costly permutation testing. This procedure thus has two
hyperparameters: the length of a search space as well as a noise threshold (in
the form of a significance level alpha) or alternatively a fixed embedding
dimensionality. In our study, we fixed these to a search space of 1 second and
a fixed embedding dimensionality of 50 (simulations) or 25 (MEG data).
Conditional MI based TE estimators that employ such multidimensional target
past operationalisations have been described to be “self prediction optimal”
and are thus referred to as TESPO (Wibral et al., 2013).

$TE_{SPO}(X,Y)_{\delta_{a}}=MI\left(X_{\delta_{a}},\left[Y_{emb}Y_{0}\right]\right)-MI\left(X_{\delta_{a}},Y_{emb}\right)$
(10)

Here, $Y_{emb}$ denotes a multidimensional embedding of the target past.

#### Whitening approaches

The suggestion we make here is to tackle the TE estimation problem with a
whitening approach. The basic idea is to use a non-uniform multidimensional
target past embedding as a model to derive a prediction about the target
present state. Analogous to (Haugh, 1976), this approach could then be used to
subtract the prediction from the observed target time series and compute
delayed MI on the residuals. This however ignores the uncertainty associated
with each prediction. We thus instead compute a time series of sample wise
conditional entropy of each sample in the observed target time series given
its past embedding. We obtain the conditional entropy of each sample $x$ given
its non-uniform embedding $x_{emb}$ by subtracting the marginal entropy of
$x_{emb}$ from the joint entropy of $x$ and $x_{emb}$:

$h\left(x\mid
x_{emb}\right)=h\left(\left[xx_{emb}\right]\right)-h\left(x_{emb}\right)$ (11)

The resulting time series is thus effectively temporally decorrelated, or
“whitened”. We define DIce between the source $X$ and the target $Y$ at an
analysis delay $\delta_{a}$ as:

$DI_{ce}(X,Y)_{\delta a}=MI\left(h\left(x\mid
x_{emb}\right)_{\delta_{a}},h\left(y\mid y_{emb}\right)\right)$ (12)

Any resulting quantifiable delayed MI between the whitened source and the
target can thus no longer stem from information that the target carries about
itself. To avoid asymmetrical temporal distortions caused by the temporal
whitening, we apply the whitening to both source and target time series.

### Partial information decomposition

In Partial information decomposition (PID, Williams & Beer, 2010), systems of
three or more variables are considered. In situations where there are two
sources and one target variable, it is the goal to quantify the “redundant”
amount of information that the two sources share about the target as well as
the “synergistic” information about the target that is only available when
considering the two sources jointly as well as the unique information about
the target that is only available from one of the sources but not the other.
From the lens of this formalism, conditional MI can be seen as the sum of two
“atoms” of PID: Unique information and synergy. TE based on conditional MI is
thus the sum of not only unique information of the source signal about the
target present but also synergy of the source signal and the target past about
the target present.

PID has its origins in co-information (McGill, 1954), which is computed as a
triple set intersection.

$\operatorname{CoI}\left(X_{\delta_{a}},Y_{emb},Y_{0}\right)=MI\left(X_{\delta_{a}},Y_{0}\right)+MI\left(Y_{emb},Y_{0}\right)-MI\left(\left[X_{\delta_{a}}Y_{emb}\right],Y_{0}\right)$
(13)

Here, negative values correspond to synergistic information, and positive
values correspond to redundant information. From the lens of PID however,
these are mere net sums of the “pure” redundancy and synergy that PID aims to
resolve. To do so, the essential first step is to define redundancy. Here, we
turned to an implementation based on “common change in surprisal” ($I_{ccs}$
Ince, 2017). This starts at the observation that local MI can be seen as the
positive or negative change in surprisal of a given value when another is
observed:

$\operatorname{mi}(x,y)=\Delta_{y}h(x)=h(x)-h(x\mid y)$ (14)

We can then also consider equally positive or negative co-information at the
local level:

$\operatorname{coi}\left(x_{\delta_{a}},y_{emb},y_{0}\right)=\operatorname{mi}\left(x_{\delta_{a}},y_{0}\right)+mi\left(y_{emb},y_{0}\right)+mi\left(\left[x_{\delta_{a}}y_{emb}\right],y_{0}\right)$
(15)

Each of such co-information terms can then either contribute to net synergy or
net redundancy. The key idea in $I_{ccs}$ is then to consider only those
positive (redundant) terms of local $coi(x_{\delta_{a}},y_{emb},y_{0})$ which
coincide with positive terms of $mi(x_{\delta_{a}},y_{0})$,
$mi(y_{emb},y_{0})$ and $mi([x_{\delta_{a}}y_{emb}],y_{0})$. Thus, only co-
information terms are counted that represent a commonly shared change in
surprisal. The resulting global redundancy can then be used to infer the other
PID atoms according to a lattice structure (Williams & Beer, 2010):

$\begin{gathered}\text{ unique
}\left(X_{\delta_{a}}\right)=MI\left(X_{\delta_{a}},Y_{0}\right)-\text{
redundancy }\\\ \text{ unique
}\left(Y_{emb}\right)=MI\left(Y_{emb},Y_{0}\right)-\text{ redundancy }\\\
\text{ synergy
}=MI\left(\left[X_{\delta_{a}}Y_{emb}\right],Y_{0}\right)-\text{ redundancy
}-\text{ unique }\left(X_{\delta_{a}}\right)-\text{ unique
}\left(Y_{emb}\right)\end{gathered}$ (16)

### Noise thresholds

To establish whether a given effect size exceeded a level that would be
expected of data not containing the effect of interest, we considered the
$95th$ percentile of distributions of effects obtained from $1000$
permutations. For this, we performed circular shifts of the source variable by
a random number of samples and then recomputed the estimators of interest. We
constrained this random amount with a minimum and a maximum, such that the
permuted data could not include instances where potential effects of the
observed data would end up within the range of delays of interest. In cases
where multiple comparisons were made, we corrected the noise threshold by
means of the family wise error rate, that is, by considering the maximum
across all conditions within a given permutation and subsequently applying the
$95th$ percentile of the resulting distribution for all conditions.

### Simulations

The basic approach to simulating narrow-band time series was to firstly sample
Gaussian white noise and subsequently apply band-pass filters to it ($3rd$
order butterworth, forward-only). Such signals $M$ were then circularly
shifted by a ground truth interaction delay. Finally, source and target time
series $X$ and $Y$ were generated by adding independent white noise of varying
amplitude to $M$. The following specific parameters were used for the
individual simulations and analyses:

Figure |  Number of samples |  Noise amplitude source |  Noise amplitude target |  Effect filter $[$Hz$]$ |  Ground truth delay $[$s$]$ |  Analysis filter $[$Hz$]$ |  Lags $[$s$]$
---|---|---|---|---|---|---|---
1 | $75000$ | $.5$ | $.25$ | $6\pm 2$ | $-.12$ | $6\pm 2$ |  $-.8$ - $+.32$
2 | $75000$ | $.5$ | $.25$ | $6\pm 2$ | $-.12$ | $6\pm 2$ |  $-.8$ - $+.32$
3 | $40000$ | $.05$ | $.025$ | $6\pm 2$ |  $-.08$ - $-.4$ | $6\pm 2$ |  $.66$ - $+.4$
4 | $40000$ |  $.001$ - $10$ |  $.001$- $10$ | $6\pm 2$ | $-.12$ | $6\pm 2$ |  $-.4$ - $+.2$
5 | $40000$ |  $0$ - $10$ |  $0$ - $10$ | $6\pm 2$ | $-.12$ | $6\pm 2$ |  $-.4$ - $+.2$
6.1 | $40000$ | $.01$ | $.01$ |  $15\pm 1$ - $13.25$ | $-.12$ |  $15\pm$1 - $13.25$ |  $-.4$ - $+.2$
6.2 | $40000$ | $.01$ | $.01$ | $20\pm 5$ | $-.12$ |  $20\pm 1$ - $15$ |  $-.4$ - $+.2$
7.1 | $40000$ | $.01$ | $.01$ |  $3$ - $50\pm 2$ | $-.12$ |  $3$ - $50\pm 2$ |  $-.4$ - $+.2$
7.2 | $40000$ | $.01$ | $.01$ |  $3$ - $47\pm 2$ | $-.12$ |  $3$ - $47\pm 2$ |  $-.4$ - $+.2$
8 | $330000$ | $na$ | $na$ | $na$ | $na$ |  $1.5$ - $16\pm 1$ |  $-1$ - $.5$
Table 1: Parameter settings in simulations. Figure 6.1 refers to panels A, B,
D and E, Figure 6.2 refers to panels C and F. Figure 7.1 refers to panels A -
D, Figure 7.2 refers to panels E - H. Effect and analysis filters are
specified in terms of the passbands as defined by centre frequencies $\pm$
bandwidths. $na$ denotes “not available”.

### MEG data and analyses

The MEG data in this study has been recorded and analysed before. For details
concerning the participants, the experimental design, the recording procedures
as well as the preprocessing, please refer to (Daube et al., 2019). In brief,
$24$ participants had listened to an audiobook of $55$ minutes duration (in 6
blocks of equal duration) while their MEG (MAGNES 3600 WH, 4D Neuroimaging,
$248$ magnetometers) had been recorded at a sampling rate of $1017.25$ Hz
(first $10$ participants) or $2034.51$ Hz (last $14$ participants). We applied
the same pre-processing steps as in the original study. This included
interpolation of artifactual channels, replacement of squid jumps with DC
patches, 4th order zero-phase high-pass filter of $.5$ Hz as well as
independent component analysis for removal of eye and heart activity (Daube et
al., 2019). Here, we then downsampled the data to a sampling rate of $100$Hz.

To estimate activity from bilateral auditory cortices (ACs), we re-used
linearly constrained minimum variance beamformer spatial filters (Van Veen et
al., 1997) as in (Daube et al., 2019). These had been optimised within a
nested cross-validation to return responses that would be maximally correlated
with linear predictions of the responses based on a combination of log-mel
spectrograms and their temporal derivative. This correlation had been
maximised with respect to position- and regularisation hyperparameters per
participant, hemisphere and fold. Here, we averaged these hyperparameters
across folds to then extract time series of activity from left and right ACs.

To define individual delta and theta frequency ranges, we considered spectra
of delayed MI between AC activity and the speech envelope. To compute these
spectra, we applied $3rd$ order butterworth forward filters to both AC
activity and the envelope. These had centre frequencies increasing from $1.5$
to $15$ Hz in steps of $.25$ Hz, and cutoff frequencies were defined as bands
of $\pm 1$ Hz width around the centre frequencies. We computed delayed MI for
a range of lags from $-1$ to $.5$s and used the maximum across lags within
each frequency. We then searched for troughs in the resulting spectra and used
the trough corresponding to the lowest frequency as the boundary between delta
and theta. In $4$ out of $48$ cases, this returned boundaries between delta
and theta that surpassed $6.5$ Hz, which we considered as unlikely a priori
(Klimesch, 1999; Wang, 2010). Therefore, in these cases we instead searched
for peaks of the spectral gradient and used the peak corresponding to the
lowest centre frequency. In practise, this corresponded to parts of the
spectra where the initial decrease plateaued. To define an upper boundary of
the theta component, we used the highest frequency at which the spectrum
surpassed the noise threshold. This resulted in theta components with a
broader bandwidth compared to the delta components. Since we had found the
bandwidth of an effect to be correlated with the effect size recovered by
DIce, we also repeated analyses of delayed MI and DIce with fixed bandwidths
for delta and theta. For this, we defined the delta component to be below
$3.5$ Hz and the theta component to be centered on a theta peak frequency,
around which we defined bandpass filters of $3$ Hz width. The theta peak
frequencies were defined as the maximum of the MI spectra within the theta
frequency band as defined previously. In $4$ out of $48$ cases, this
corresponded to the upper boundary of the delta band. In these cases, we
defined the theta centre frequency to be $1.5$ Hz above the delta-theta
boundary.

To extract time series corresponding to the passbands as defined above, we
used $3rd$ order forward butterworth filters, which we applied to both the
envelope and the AC activity. We then subjected the resulting time series to
the computation of delayed MI and DIce. We statistically evaluated the nth
observation of recovered log-transformed effect sizes and the recovered delays
$r$ at the maxima across lags of each participant $p$, hemisphere $h$,
frequency band $b$ and measure $m$ by fitting Bayesian linear distributional
models as implemented in the brms package (Bürkner, 2017). These models can be
summarised with the following formula:

$\displaystyle r_{[n]}$ $\displaystyle\sim
N\left(\eta_{\mu[n]},\exp\left(\eta_{\sigma[n]}\right)\right)$ (17)
$\displaystyle\eta_{\mu[n]}$
$\displaystyle\sim\beta_{\mu_{p[n]}}+\beta_{\mu_{p:h[n]}}+\beta_{\mu_{m:b[n]}}$
$\displaystyle\left(\beta_{\mu_{p[n]}},\beta_{\mu_{p:h[n]^{\prime}}}\beta_{\mu_{m:b[n]}}\right)$
$\displaystyle\sim N(0,v)$ $\displaystyle\eta_{\sigma[n]}$
$\displaystyle\sim\beta_{\left.\sigma_{h[n]}\right]}+\beta_{\sigma_{m:b[n]}}$
$\displaystyle\left(\beta_{\left.\sigma_{h[n]}\right]},\beta_{\sigma_{m:b[n]}}\right)$
$\displaystyle\sim N(0,10)$

## Supplementary Material

Figure S1: Spectra of MI for left and right auditory cortices for each
individual participant (related to figure 9). Individual boundaries of delta
and theta bands are overlaid as filled circles, theta peak frequencies (used
for figure S3) are overlaid as empty circles.  Figure S2: Caption on following
page. Figure S3: Results obtained from source level MEG recordings with
constrained frequencies (related to figure 9). Same as figure 9, but computed
with constrained bandwidths of analysis filters. Delta band was fixed by means
of a low-pass filter common across all participants with a cutoff frequency of
$3.5$Hz. Theta frequency was defined as a $3$Hz wide band centered on
individual theta centre frequencies (see figure S1). The increased DI/MI ratio
found with individualised band definitions (where often theta was defined as a
wider band than delta) was robust to the constrained analysis bandwidths,
making it unlikely that bandwidth alone can fully account for the difference
in DI/MI ratio between delta and theta bands. Also see figure S4.  Figure S4:
Comparisons of posterior distributions of main effects of combinations of
frequency bands and dependency measures from Bayesian linear modeling of the
raw dependency and lag estimates (related to figure 9). In each cell in the
matrices, the greyscale colour denotes the fraction of samples of the
combination of frequency band and dependency measure referenced on the y-axis
that is larger than the combination referenced on the x-axis (testing a
hypothesis). Top row reports results corresponding to figure 9, bottom row
reports results corresponding to figure S3.

## Contributions

C.D., J.G. and R.A.A.I. conceived of and designed the experiment.

C.D. collected and analyzed the data.

C.D. and R.A.A.I. contributed analytic tools.

C.D. wrote the manuscript.

C.D., J.G. and R.A.A.I. edited the manuscript.

## Acknowledgements

R.A.A.I. was supported by the Wellcome Trust [214120/Z/18/Z].

## References

  * Ahissar et al. (2001) Ahissar, E., Nagarajan, S., Ahissar, M., Protopapas, A., Mahncke, H., & Merzenich, M. M. (2001).  Speech comprehension is correlated with temporal response patterns recorded from auditory cortex.  Proceedings of the National Academy of Sciences, 98(23):13367–13372.
  * Barnett et al. (2009) Barnett, L., Barrett, A. B., & Seth, A. K. (2009).  Granger causality and transfer entropy are equivalent for Gaussian variables.  Physical review letters, 103(23):238701.
  * Barnett & Seth (2011) Barnett, L. & Seth, A. K. (2011).  Behaviour of Granger causality under filtering: theoretical invariance and practical application.  Journal of neuroscience methods, 201(2):404–419.
  * Barrett (2015) Barrett, A. B. (2015).  Exploration of synergistic and redundant information sharing in static and dynamical Gaussian systems.  Physical Review E, 91(5):052802.
  * Bassett & Bullmore (2006) Bassett, D. S. & Bullmore, E. (2006).  Small-world brain networks.  The neuroscientist, 12(6):512–523.
  * Bastos & Schoffelen (2016) Bastos, A. M. & Schoffelen, J.-M. (2016).  A tutorial review of functional connectivity analysis methods and their interpretational pitfalls.  Frontiers in systems neuroscience, 9:175.
  * Besserve et al. (2010) Besserve, M., Schölkopf, B., Logothetis, N. K., & Panzeri, S. (2010).  Causal relationships between frequency bands of extracellular signals in visual cortex revealed by an information theoretic analysis.  Journal of computational neuroscience, 29(3):547–566.
  * Brodbeck et al. (2018) Brodbeck, C., Hong, L. E., & Simon, J. Z. (2018).  Rapid Transformation from Auditory to Linguistic Representations of Continuous Speech.  Current Biology, 28:3976–3983.
  * Brookes et al. (2011) Brookes, M. J., Woolrich, M., Luckhoo, H., Price, D., Hale, J. R., Stephenson, M. C., Barnes, G. R., Smith, S. M., & Morris, P. G. (2011).  Investigating the electrophysiological basis of resting state networks using magnetoencephalography.  Proceedings of the National Academy of Sciences, 108(40):16783–16788.
  * Bürkner (2017) Bürkner, P. C. (2017).  brms: An R Package for Bayesian Multilevel Models Using Stan.  Journal of Statistical Software, 80(1).
  * Caucheteux et al. (2021) Caucheteux, C., Gramfort, A., & King, J.-R. (2021).  Long-range and hierarchical language predictions in brains and algorithms.  arXiv preprint arXiv:2111.14232.
  * Chung et al. (2020) Chung, Y.-A., Tang, H., & Glass, J. (2020).  Vector-quantized autoregressive predictive coding.  arXiv preprint arXiv:2005.08392.
  * Cliff et al. (2021) Cliff, O. M., Novelli, L., Fulcher, B. D., Shine, J. M., & Lizier, J. T. (2021).  Assessing the significance of directed and multivariate measures of linear dependence between time series.  Physical Review Research, 3(1):013145.
  * Crosse et al. (2016) Crosse, M. J., DiLiberto, G. M., Bednar, A., & Lalor, E. C. (2016).  The Multivariate Temporal Response Function (mTRF) Toolbox: A MATLAB Toolbox for Relating Neural Signals to Continuous Stimuli.  frontiers in Human Neuroscience, 10(604).
  * Daube et al. (2019) Daube, C., Ince, R. A. A., & Gross, J. (2019).  Simple Acoustic Features Can Explain Phoneme-Based Predictions of Cortical Responses to Speech.  Current Biology, 29(12):1924–1937.e9.  doi: 10.1016/j.cub.2019.04.067.  URL https://www.sciencedirect.com/science/article/pii/S0960982219304968.
  * de-Wit et al. (2016) de-Wit, L., Alexander, D., Ekroll, V., & Wagemans, J. (2016).  Is neuroimaging measuring information in the brain?  Psychonomic Bulletin & Reviews, 23:1415—1428.
  * Di Liberto et al. (2015) Di Liberto, G. M., O’Sullivan, J. A., & Lalor, E. C. (2015).  Low-Frequency Cortical Entrainment to Speech Reflects Phoneme-Level Processing.  Current Biology, 25:2457–2465.
  * Ding et al. (2016) Ding, N., Melloni, L., Zhang, H., Tian, X., & Poeppel, D. (2016).  Cortical tracking of hierarchical linguistic structures in connected speech.  Nature Neuroscience, 19(1):158–164.
  * Ding & Simon (2012) Ding, N. & Simon, J. Z. (2012).  Emergence of neural encoding of auditory objects while listening to competing speakers.  Proceedings of the National Academy of Sciences of the United States of America, 109(29):11854–11859.
  * Ding & Simon (2014) Ding, N. & Simon, J. Z. (2014).  Cortical entrainment to continuous speech: functional roles and interpretations.  Frontiers in human neuroscience, 8:311.
  * Donhauser & Baillet (2020) Donhauser, P. W. & Baillet, S. (2020).  Two distinct neural timescales for predictive speech processing.  Neuron, 105(2):385–393.
  * Donoghue et al. (2020) Donoghue, T., Haller, M., Peterson, E. J., Varma, P., Sebastian, P., Gao, R., Noto, T., Lara, A. H., Wallis, J. D., Knight, R. T., et al. (2020).  Parameterizing neural power spectra into periodic and aperiodic components.  Nature neuroscience, 23(12):1655–1665.
  * Etard & Reichenbach (2019) Etard, O. & Reichenbach, T. (2019).  Neural speech tracking in the theta and in the delta frequency band differentially encode clarity and comprehension of speech in noise.  Journal of Neuroscience, 39(29):5750–5759.
  * Faes et al. (2011) Faes, L., Nollo, G., & Porta, A. (2011).  Information-based detection of nonlinear Granger causality in multivariate processes via a nonuniform embedding technique.  Physical Review E, 83(5):051112.
  * Faes et al. (2017) Faes, L., Nollo, G., Stramaglia, S., & Marinazzo, D. (2017).  Multiscale granger causality.  Physical Review E, 96(4):042150.
  * Florin et al. (2010) Florin, E., Gross, J., Pfeifer, J., Fink, G. R., & Timmermann, L. (2010).  The effect of filtering on Granger causality based multivariate causality measures.  Neuroimage, 50(2):577–588.
  * Fries (2015) Fries, P. (2015).  Rhythms for cognition: communication through coherence.  Neuron, 88(1):220–235.
  * Friston (2005) Friston, K. (2005).  A theory of cortical responses.  Philosophical transactions of the Royal Society B: Biological sciences, 360(1456):815–836.
  * Gerster et al. (2021) Gerster, M., Waterstraat, G., Litvak, V., Lehnertz, K., Schnitzler, A., Florin, E., Curio, G., & Nikulin, V. (2021).  Separating neural oscillations from aperiodic 1/f activity: challenges and recommendations.  bioRxiv.
  * Giordano et al. (2016) Giordano, B. L., Ince, R. A. A., Gross, J., Schyns, P. G., Panzeri, S., & Kayser, C. (2016).  Contributions of local speech encoding and functional connectivity to audio-visual speech perception.  eLife, 6:e24763(DOI: 10.7554/eLife.24763).
  * Giraud & Poeppel (2012) Giraud, A. L. & Poeppel, D. (2012).  Cortical oscillations and speech processing: emerging computational principles and operations.  Nature Neuroscience, 15:511–517.
  * Granger (1969) Granger, C. W. (1969).  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica: journal of the Econometric Society, pages 424–438.
  * Gross et al. (2013) Gross, J., Hoogenboom, N., Thut, G., Schyns, P. G., Panzeri, S., Belin, P., & Garrod, S. (2013).  Speech Rhythms and Multiplexed Oscillatory Sensory Coding in the Human Brain.  PLoS Biology, 11(12).
  * Gross et al. (2021) Gross, J., Kluger, D. S., Abbasi, O., Chalas, N., Steingräber, N., Daube, C., & Schoffelen, J.-M. (2021).  Comparison of undirected frequency-domain connectivity measures for cerebro-peripheral analysis.  NeuroImage, 245:118660.
  * Haugh (1976) Haugh, L. D. (1976).  Checking the independence of two covariance-stationary time series: a univariate residual cross-correlation approach.  Journal of the American Statistical Association, 71(354):378–385.
  * Heilbron et al. (2021) Heilbron, M., Armeni, K., Schoffelen, J.-M., Hagoort, P., & de Lange, F. P. (2021).  A hierarchy of linguistic predictions during natural language comprehension.  bioRxiv, pages 2020–12.
  * Henry & Obleser (2012) Henry, M. J. & Obleser, J. (2012).  Frequency modulation entrains slow neural oscillations and optimizes human listening behavior.  Proceedings of the National Academy of Sciences, 109(49):20095–20100.
  * Hertrich et al. (2012) Hertrich, I., Dietrich, S., Trouvain, J., Moos, A., & Ackermann, H. (2012).  Magnetic brain activity phase-locked to the envelope, the syllable onsets, and the fundamental frequency of a perceived speech signal.  Psychophysiology, 49:322–334.
  * Hyafil et al. (2015) Hyafil, A., Fontolan, L., Kabdebon, C., Boris, G., & Giraud, A. L. (2015).  Speech encoding by coupled cortical theta and gamma oscillations.  eLife, 4(e06213).
  * Ince (2017) Ince, R. (2017).  Measuring Multivariate Redundant Information with Pointwise Common Change in Surprisal.  Entropy, 19(7):318.  doi: 10.3390/e19070318.  URL http://www.mdpi.com/1099-4300/19/7/318.
  * Ince et al. (2015) Ince, R. A., Van Rijsbergen, N. J., Thut, G., Rousselet, G. A., Gross, J., Panzeri, S., & Schyns, P. G. (2015).  Tracing the flow of perceptual features in an algorithmic brain network.  Scientific reports, 5(1):1–17.
  * Ince et al. (2017) Ince, R. A. A., Giordano, B. L., Kayser, C., Rousselet, G. A., Gross, J., & Schyns, P. G. (2017).  A Statistical Framework for Neuroimaging Data Analysis Based on Mutual Information Estimated via a Gaussian Copula.  Human Brain Mapping, 38(3):1541–1573.
  * Ince et al. (2021) Ince, R. A. A., Kay, J. W., & Schyns, P. G. (2021).  Bayesian inference of population prevalence.  eLife, 10(e62461).  doi: 10.7554/eLife.62461.  URL https://elifesciences.org/articles/62461.
  * Jain et al. (2021) Jain, S., Mahto, S., Turek, J. S., Vo, V. A., LeBel, A., & Huth, A. G. (2021).  Interpretable multi-timescale models for predicting fMRI responses to continuous natural speech.  bioRxiv, pages 2020–10.
  * James et al. (2016) James, R. G., Barnett, N., & Crutchfield, J. P. (2016).  Information flows? A critique of transfer entropies.  Physical review letters, 116(23):238701.
  * Jones & Boltz (1989) Jones, M. R. & Boltz, M. (1989).  Dynamic attending and responses to time.  Psychological review, 96(3):459.
  * Kayser et al. (2015) Kayser, C., Wilson, C., Safaai, H., Sakata, S., & Panzeri, S. (2015).  Rhythmic auditory cortex activity at multiple timescales shapes stimulus–response gain and background firing.  Journal of Neuroscience, 35(20):7750–7762.
  * Klimesch (1999) Klimesch, W. (1999).  EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis.  Brain research reviews, 29(2-3):169–195.
  * Koskinen et al. (2020) Koskinen, M., Kurimo, M., Gross, J., Hyvärinen, A., & Hari, R. (2020).  Brain activity reflects the predictability of word sequences in listened continuous speech.  NeuroImage, 219:116936.
  * Kulasingham et al. (2020) Kulasingham, J. P., Brodbeck, C., Presacco, A., Kuchinsky, S. E., Anderson, S., & Simon, J. Z. (2020).  High gamma cortical processing of continuous speech in younger and older listeners.  Neuroimage, 222:117291.
  * Lakatos et al. (2019) Lakatos, P., Gross, J., & Thut, G. (2019).  A new unifying account of the roles of neuronal entrainment.  Current Biology, 29(18):R890–R905.
  * Lakatos et al. (2008) Lakatos, P., Karmos, G., Mehta, A. D., Ulbert, I., & Schroeder, C. E. (2008).  Entrainment of neuronal oscillations as a mechanism of attentional selection.  science, 320(5872):110–113.
  * Lakhotia et al. (2021) Lakhotia, K., Kharitonov, E., Hsu, W.-N., Adi, Y., Polyak, A., Bolte, B., Nguyen, T.-A., Copet, J., Baevski, A., Mohamed, A., et al. (2021).  Generative spoken language modeling from raw audio.  arXiv preprint arXiv:2102.01192.
  * Lalor et al. (2009) Lalor, E. C., Power, A. J., Reilly, R. B., & Foxe, J. J. (2009).  Resolving precise temporal processing properties of the auditory system using continuous stimuli.  Journal of neurophysiology, 102(1):349–359.
  * Lobier et al. (2014) Lobier, M., Siebenhühner, F., Palva, S., & Palva, J. M. (2014).  Phase transfer entropy: a novel phase-based measure for directed connectivity in networks coupled by oscillatory interactions.  Neuroimage, 85:853–872.
  * Makeig et al. (2002) Makeig, S., Westerfield, M., Jung, T.-P., Enghoff, S., Townsend, J., Courchesne, E., & Sejnowski, T. J. (2002).  Dynamic brain sources of visual evoked responses.  Science, 295(5555):690–694.
  * McGill (1954) McGill, W. J. (1954).  Multivariate information transmission.  Psychometrika, 19:97–116.
  * Mehler & Kording (2018) Mehler, D. M. A. & Kording, K. P. (2018).  The lure of causal statements: Rampant mis-inference of causality in estimated connectivity.  arXiv e-prints, pages arXiv–1812.
  * Mell et al. (2021) Mell, M. M., St-Yves, G., & Naselaris, T. (2021).  Voxel-to-voxel predictive models reveal unexpected structure in unexplained variance.  NeuroImage, page 118266.
  * Michalareas et al. (2016) Michalareas, G., Vezoli, J., Van Pelt, S., Schoffelen, J.-M., Kennedy, H., & Fries, P. (2016).  Alpha-beta and gamma rhythms subserve feedback and feedforward influences among human visual cortical areas.  Neuron, 89(2):384–397.
  * Młynarski & Hermundstad (2018) Młynarski, W. F. & Hermundstad, A. M. (2018).  Adaptive coding for dynamic sensory inference.  Elife, 7:e32055.
  * Morillon & Baillet (2017) Morillon, B. & Baillet, S. (2017).  Motor origin of temporal predictions in auditory attention.  Proceedings of the National Academy of Sciences, 114(42):E8913–E8921.
  * Naselaris et al. (2011) Naselaris, T., Kay, K. N., Nishimoto, S., & Gallant, J. L. (2011).  Encoding and decoding in fMRI.  NeuroImage, 56(2):400–410.  doi: 10.1016/j.neuroimage.2010.07.073.  URL http://www.sciencedirect.com/science/article/pii/S1053811910010657.
  * Obleser & Kayser (2019) Obleser, J. & Kayser, C. (2019).  Neural entrainment and attentional selection in the listening brain.  Trends in cognitive sciences, 23(11):913–926.
  * O’Sullivan et al. (2015) O’Sullivan, J. A., Power, A. J., Mesgarani, N., Rajaram, S., Foxe, J. J., Shinn-Cunningham, B. G., Slaney, M., Shamma, S. A., & Lalor, E. C. (2015).  Attentional selection in a cocktail party environment can be decoded from single-trial EEG.  Cerebral cortex, 25(7):1697–1706.
  * Park et al. (2015) Park, H., Ince, R. A., Schyns, P. G., Thut, G., & Gross, J. (2015).  Frontal top-down signals increase coupling of auditory low-frequency oscillations to continuous speech in human listeners.  Current Biology, 25(12):1649–1653.
  * Pinzuti et al. (2020) Pinzuti, E., Wollstadt, P., Gutknecht, A., Tüscher, O., & Wibral, M. (2020).  Measuring spectrally-resolved information transfer.  PLoS Computational Biology, 16(12):e1008526.
  * Ragwitz & Kantz (2002) Ragwitz, M. & Kantz, H. (2002).  Markov models from data by simple nonlinear time series predictors in delay embedding spaces.  Physical Review E, 65(5):056201.
  * Rao & Ballard (1999) Rao, R. P. & Ballard, D. H. (1999).  Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.  Nature neuroscience, 2(1):79–87.
  * Ritter et al. (1968) Ritter, W., Vaughan Jr, H. G., & Costa, L. D. (1968).  Orienting and habituation to auditory stimuli: a study of short terms changes in average evoked responses.  Electroencephalography and clinical Neurophysiology, 25(6):550–556.
  * Sayers et al. (1974) Sayers, B. M., Beagley, H. A., & Henshall, W. R. (1974).  The mechanism of auditory evoked EEG responses.  Nature, 247(5441):481–483.
  * Schmitt et al. (2021) Schmitt, L.-M., Erb, J., Tune, S., Rysop, A., Hartwigsen, G., & Obleser, J. (2021).  Predicting speech from a cortical hierarchy of event-based time scales.  Science Advances, 7(49).
  * Schnitzler & Gross (2005) Schnitzler, A. & Gross, J. (2005).  Normal and pathological oscillatory communication in the brain.  Nature reviews neuroscience, 6(4):285–296.
  * Schoffelen et al. (2017) Schoffelen, J.-M., Hultén, A., Lam, N., Marquand, A. F., Uddén, J., & Hagoort, P. (2017).  Frequency-specific directed interactions in the human brain network for language.  Proceedings of the National Academy of Sciences, 114(30):8083–8088.
  * Schreiber (2000) Schreiber, T. (2000).  Measuring information transfer.  Physical review letters, 85(2):461.
  * Schroeder & Lakatos (2009) Schroeder, C. E. & Lakatos, P. (2009).  Low-frequency neuronal oscillations as instruments of sensory selection.  Trends in neurosciences, 32(1):9–18.
  * Takens (1981) Takens, F. (1981).  Detecting strange attractors in turbulence.  In Dynamical systems and turbulence, Warwick 1980, pages 366–381. Springer.
  * Van Veen et al. (1997) Van Veen, B. D., van Drongelen, W., Yuchtman, M., & Suzuki, A. (1997).  Localization of Brain Electrical Activity via Linearly Constrained Minimum Variance Spatial Filtering.  IEEE Transactions on Biomedical Engineering, 44(9):867–880.
  * Vlachos & Kugiumtzis (2010) Vlachos, I. & Kugiumtzis, D. (2010).  Nonuniform state-space reconstruction and coupling detection.  Physical Review E, 82(1):016207.
  * Wang (2010) Wang, X.-J. (2010).  Neurophysiological and computational principles of cortical rhythms in cognition.  Physiological reviews, 90(3):1195–1268.
  * Wibral et al. (2013) Wibral, M., Pampu, N., Priesemann, V., Siebenhühner, F., Seiwert, H., Lindner, M., Lizier, J. T., & Vicente, R. (2013).  Measuring information-transfer delays.  PloS one, 8(2):e55809.
  * Williams & Beer (2010) Williams, P. L. & Beer, R. D. (2010).  Nonnegative Decomposition of Multivariate Information.  arXiv, 1004.2515.
  * Wollstadt et al. (2017) Wollstadt, P., Sellers, K. K., Rudelt, L., Priesemann, V., Hutt, A., Fröhlich, F., & Wibral, M. (2017).  Breakdown of local information processing may underlie isoflurane anesthesia effects.  PLoS computational biology, 13(6):e1005511.
  * Wöstmann et al. (2017) Wöstmann, M., Fiedler, L., & Obleser, J. (2017).  Tracking the signal, cracking the code: Speech and speech comprehension in non-invasive human electrophysiology.  Language, Cognition and Neuroscience, 32(7):855–869.
  * Zan et al. (2020) Zan, P., Presacco, A., Anderson, S., & Simon, J. Z. (2020).  Exaggerated cortical representation of speech in older listeners: mutual information analysis.  Journal of Neurophysiology, 124(4):1152–1164.