Efficient transfer entropy analysis of non-stationary neural time series

Introduction

• Information transfer is not equivalent to causal interaction, and distinguishing transfer, storage, and modification can help identify where and when neural information processing occurs.The passages note that causal interactions can serve all three information-processing components, whereas TE has become a widely used measure of information transfer [19– 41].Source: Introduction · pp. 1–2
• Multivariate TE can condition on other system processes, but limited samples hinder naive estimation; proposed approximations include iterative inclusion methods and delay reconstruction with graphical pruning.The authors describe the graphical method as potentially better suited to limited data, while the other approaches may come closer to ground truth; they consider the multivariate obstacle practically addressed when enough data are available.Source: Introduction · p. 2
• The stationarity requirement is avoidable when multiple realizations are available; it is needed to estimate probability density functions from a single realization when repetitions are unavailable.Neuroscience experiments often obtain realizations by repeating trials under as-constant-as-possible conditions.Source: Introduction · p. 2
• Ensemble-based TE estimation requires many surrogate datasets, each with a TE estimate, and the MEG analysis described would otherwise take 8200 h for 15 subjects and one condition.The nearest-neighbor KSG estimator is described as precise but computationally intensive, making a substantial speed-up necessary for practical use.Source: Introduction · p. 2
• Repeated experimental trials provide an ensemble for time-resolved TE estimation without assuming stationarity, but estimation must pool across trials and use many surrogate datasets for statistical testing.Gomez-Herrero and colleagues [55] demonstrated the ensemble approach theoretically; unlike stationary-data analysis, it does not estimate TE per trial or reuse one set of trial-based surrogates.Source: Introduction · p. 2
• A GPU implementation parallelizing neighbor searches reduced computation time 50-fold versus the previously used CPU-optimized TE search, making ensemble-based analysis feasible for typical neural datasets.The implementation parallelized the computationally heaviest part of binless TE estimation and enabled time-resolved TE estimation from non-stationary neural time series.Source: Introduction · p. 2

Background

• TE estimation maps observed time-series data to an information-theoretic quantity by estimating the joint and conditional probability distributions of random-process variables.The background formalizes a recorded time series as a realization of a random process and calls the data-to-quantity mapping an estimator.Source: Background · pp. 2–3

Stationarity and non-stationarity in experimental time series

• Stationary processes allow observations to be pooled over time, but changing probability distributions in non-stationary processes make temporal pooling inapplicable.Repeated observations aligned by time relative to repetition onset can instead be pooled across repetitions when their distributions agree.Source: Stationarity and non-stationarity in… · p. 3
• Because physical copies are typically unavailable in neuroscience, repeated stimulus presentations provide trials for ensemble estimation; Gomez-Herrero and colleagues proposed this approach for cyclostationary processes [55].The paper treats experimental trials as repetitions of a repeatable process.Source: Stationarity and non-stationarity in… · p. 3

Transfer entropy estimation from an ensemble of time series

• Ensemble TE extends prior functionals to repeatable or cyclostationary non-stationary data, avoids temporal pooling, and enables time-resolved estimates with self-prediction optimality [55,19,20,53].The functional explicitly obtains realizations through ensemble evaluation over repetitions.Source: Transfer entropy estimation from an… · pp. 3–4
• The self-prediction-optimality functional uses the source past at an assumed transfer delay u alongside the target’s present and past, and contains Schreiber’s TE as a special case at u=1 [14,53].The functional is an extension of earlier TE formulations [19,20,53].Source: Transfer entropy estimation from an… · p. 3
• Practical estimation reconstructs finite delay embeddings that approximate the processes’ states, with embedding parameters chosen to optimize a local predictor and avoid TE overestimation [53,71].The text notes that finite-order Markov approximations can work for practical purposes even for high-dimensional stochastic systems.Source: Transfer entropy estimation from an… · p. 5
• The modified KSG estimator fixes k in the highest-dimensional space and projects neighbor-search distances into lower-dimensional spaces, searching across repetitions [56,57,72,73,74].This avoids the differing spatial scales that would result from separately fixing neighbor counts in each space.Source: Transfer entropy estimation from an… · p. 5
• When repetitions alone are insufficient, combining ensemble evaluation with time-window pooling provides more observations, while window width controls temporal resolution.The method counts neighbors within a time window around each target time.Source: Transfer entropy estimation from an… · p. 5
• Because nearest-neighbor searches are costly, the implementation parallelizes them on GPUs and tests original TE against surrogates generated by shuffling target repetitions.The implementation uses NVIDIA CUDA code to run thousands of searches in parallel; TRENTOOL’s workflow passes original and surrogate data to the GPU before statistical testing.Source: Transfer entropy estimation from an… · p. 6, Figure 2 · p. 6

Main analysis work flow in TRENTOOL

• The practical TRENTOOL 3.0 workflow has two stages: data preparation optimizes embedding parameters, then ensemble-based TE estimation is performed for source and target processes.For each process pair, X is the source and Y is the target.Source: Main analysis work flow in TRENTOOL · pp. 6–7
• TE estimation embeds and concatenates original and shuffled-target repetitions, searches the resulting data chunks in parallel on the GPU, and calculates TE from neighbor counts using the KSG estimator [57].The original and surrogate search spaces are independent, enabling their computations to run in parallel.Source: Main analysis work flow in TRENTOOL · p. 7

TE calculation and statistical testing against surrogate data

• Original TE is tested against shuffled-target surrogate TE values; the p-value is the proportion of surrogate values equal to or larger than the original, compared with critical alpha [56,76].Surrogates are formed by randomly permuting repetitions of the target process.Source: Main analysis work flow in TRENTOOL · p. 7

Reconstruction of information transfer delays

• Transfer delay is reconstructed by estimating TE across candidate u values and selecting the value that maximizes TE; delay reconstruction also tested implementation correctness and robustness.The delay between X and Y is denoted d_XY.Source: Main analysis work flow in TRENTOOL · p. 7

Parallelized nearest neighbor searches

• The implementation handles multiple problem instances simultaneously, addressing TE workloads with numerous small instances and supporting variable data configurations.The authors designed it to process chunks of embedded original and surrogate data in parallel and to make better use of GPU resources than handling one instance at a time.Source: Implementation of the GPU algorithm · p. 8

Input data

• For TE estimation, one GPU call searches chunks for the original time-series pair and its surrogate datasets, using temporal pooling, ensemble pooling, or both.The chunks are search spaces formed from embedded data points for a given pair of time series and chosen u and t values.Source: Implementation of the GPU algorithm · p. 8

Core algorithm

• Each CUDA thread brute-force searches one reference point within its chunk, while searches within and across chunks run in parallel.Figure 4 depicts one thread block per chunk, with data prepared on the CPU and stored in GPU global memory as Structures of Arrays.Source: Figure 4 · p. 8, Implementation of the GPU algorithm · p. 9

Low-level implementation details

• To improve GPU utilization and reduce memory latency, the implementation maximizes occupancy, coalesces memory accesses, stores points as Structures of Arrays, and tracks neighbors in shared memory.The implementation uses blocks of 512 threads.Source: Implementation of the GPU algorithm · p. 9

Performance

• GPU execution was faster than CPU execution by factors 22, 33 and 50, respectively, across three GPU devices of varying computing power.Figure 5 reports combined execution times for serial and parallel k-nearest-neighbor and range searches as input size varies.Source: Figure 5 · p. 9

Evaluation

• The evaluation tests whether the algorithm’s speedup is sufficient for real-world neural datasets.This assesses whether the method can be applied to such datasets.Source: Evaluation · p. 9
• The evaluation checks correctness on simulated data with known ground truth.Source: Evaluation · p. 9
• The evaluation tests the algorithm’s robustness when sample sizes are limited.Source: Evaluation · p. 9
• The evaluation assesses whether the algorithm produces plausible results on a neural example dataset.Source: Evaluation · p. 9

Ethics statement

• The dataset was taken from an experiment described in [86]; all subjects gave written informed consent, and the local ethics committee at Johann Wolfgang Goethe University approved the study.Source: Ethics statement · p. 9

• The benchmark compared GPU kNNS and range searches with serial TSTOOL implementations, while holding preprocessing constant to isolate search performance.Input sizes increased by duplicating chunks from simulated coupled Lorenz-system data; the GPU searched chunks in parallel, while the CPU processed them iteratively.Source: Evaluation of computational speedup · pp. 9–10
• GPU searches were faster than CPU searches by factors of 22 on the Tesla C2075, 33 on the GTX 580, and 50 on the GTX Titan (Figure 5).The CPU execution times were extrapolated to larger chunk counts and compared with measured GPU times.Source: Evaluation of computational speedup · p. 10
• Across 200 CPU runs, mean TSTOOL execution times were 1.26 s for kNNS and 24.1 s for range search.The serial timings were extrapolated to higher chunk counts; GPU tests measured the proposed searches on three NVIDIA devices.Source: Evaluation of computational speedup · p. 10
• A representative workload of about 6 million chunks was estimated to take 4.8 weeks on one GTX Titan GPU versus 240 weeks on one CPU.The estimate uses 24.1 s/50 for a typical test dataset; scanning information-transfer delays may require parallelization over multiple GPUs.Source: Evaluation of computational speedup · p. 10

Evaluation on Lorenz systems

• The Lorenz test varied unidirectional X→Y coupling over time while holding other parameters constant, with coupling limited to approximately 1000–2000 ms and a simulated 45 ms delay.The analysis used 150 repetitions of 3000 data points and evaluated three 250 ms windows.Source: Evaluation on Lorenz systems · pp. 10–11
• For Lorenz systems, significant transfer entropy appeared only during the coupling window, with a 49 ms reconstructed delay versus the true 45 ms delay (error under 10%).No significant transfer was found during the uncoupled intervals.Source: Figure 6 · p. 10, Evaluation on Lorenz systems · p. 11

Evaluation on autoregressive processes

• The AR(1) evaluation tested three coupling patterns—unidirectional onset, a two-step strength increase, and bidirectional onset—using hyperbolic-tangent modulation for smooth transitions.Each scenario used 50 trials of 3000 ms and 300 ms analysis windows; scenarios (2) and (3) were evaluated in eight windows.Source: Evaluation on autoregressive processes · pp. 11, 13
• In AR(1) scenario (1), significant transfer entropy emerged in windows 3 and 4, with maximum values at 7 ms and 9 ms, respectively.The first two windows showed no significant transfer; the delay reconstruction was less accurate during the weaker 0.8–1.1 s transition, possibly because data were limited.Source: Evaluation on autoregressive processes · p. 13
• In AR(1) scenario (2), transfer entropy was absent before coupling onset and detected in windows 4–8; reconstructed delays were 10, 4, and 7 ms in windows 4–6, then 10 ms in windows 7–8.A false positive occurred in the Y→X direction in window 6.Source: Figure 7 · p. 12, Evaluation on autoregressive processes · p. 13
• In AR(1) scenario (3), X→Y transfer was detected after its onset, with the correct 10 ms delay in windows 6–8; Y→X transfer was detected after its onset with the correct 20 ms delay in windows 7–8.Neither direction showed detected transfer before its coupling onset.Source: Figure 7 · p. 12, Evaluation on autoregressive processes · p. 15

Evaluation of the robustness of ensemble-based TE-estimation

• 10,000 and 30,000 data points reconstructed the simulated 45 ms transfer delay within 5%, while 5,000, 2,000, and 500 points yielded delays of 48, 50, and 54 ms.The test used up to 300 repetitions of 300 ms at 1000 Hz and assessed TE across assumed delays from 30 to 60 ms.Source: Evaluation of the robustness of… · p. 15

Evaluation on neural time series from magnetoencephalography

• The evaluation used MEG data from 15 healthy subjects performing a face-detection task, retaining repetitions with correctly identified face conditions.Subjects viewed degraded face or scrambled stimuli and indicated face detection by button press.Source: Evaluation on neural time series from… · p. 15
• The source analysis used a frequency-domain beamformer at 80 Hz to compare face and non-face conditions, after excluding artifact-contaminated epochs and baseline-correcting the data.The beamformer used common filters and 5% regularization for the short analysis windows.Source: Evaluation on neural time series from… · p. 15
• 5.00 is listed in all three displayed entries for cITG→FFA and MiFG→Premotor in the MEG delay table.The table also includes region pairs with missing entries and other listed delays, including 5.83 for STG→Premotor.Source: Table 2 · p. 16, Evaluation on neural time series from… · p. 16

Results

• Time-resolved GPU-based TE analysis found significant group-level information transfer changing over time; mean reconstructed delays for significant transfers were 5.00 ms (Premotor STG), 5.00 ms in all three windows (Premotor MiFG), 5.25 ms (Cing STG), and 5.67 ms (Cing FFA).Source: V1 SMG - - 5.00 · p. 16
• Fourteen sources with differential spectral power between conditions at p<0.01 were selected across occipital, parietal, temporal, and frontal cortices for broadband 10–150-Hz beamformer reconstruction.Exact anatomical locations are referenced in [86].Source: V1 SMG - - 5.00 · p. 16
• TE was estimated from beamformer source time courses by mixed pooling over repetitions and embedded time points in three non-overlapping 150-ms windows.The windows were 0–150 ms, 150–300 ms, and 300–450 ms.Source: V1 SMG - - 5.00 · p. 16
• Significant-transfer delays were scanned from 5 to 17 ms in 2-ms steps, then corrected for cascade and common-drive effects using a graph-based post-hoc method.The delay scan followed the approach in [53], and the correction method was proposed in [54].Source: V1 SMG - - 5.00 · p. 16

Efficient transfer entropy estimation from an ensemble of time series

• GPU parallelization makes the ensemble method proposed by [55] viable for short-timescale neural information-transfer analysis, implemented in the open-source TRENTOOL 3.0 toolbox.Source: Efficient transfer entropy estimation… · p. 17
• The implementation may also support large datasets in physiology [42–44], climatology [94,95], finance [45,96], and cellular automata [48], beyond neuroscience.Source: Efficient transfer entropy estimation… · p. 17

Notes on the practical application of the ensemble method for TE estimation

• Simulated data yielded reliable information-transfer reconstruction despite non-stationarity, while MEG analysis was practical and time-resolved findings aligned with prior work [86,92,93].Source: Notes on the practical application of… · p. 17
• One MEG subject and analysis window took 93 hours on average; serially analyzing 15 subjects across three windows would take about six months, reducible with multiple GPUs.Source: Notes on the practical application of… · p. 17
• The MEG example yielded 8,800 embedded data points per 150-ms window, whereas simulated TE reconstruction was correct from 10,000 points; shorter windows require more repetitions.Data available for estimation depends on window size, repetitions, and embedding parameters; the authors advise checking post-embedding data volume when designing experiments.Source: Notes on the practical application of… · p. 17
• Appropriate window and embedding choices are crucial: oversized embedding dimensions waste data and may weaken detection in noisy data, so an optimal dimension such as by Ragwitz’ criterion [71] is important [56].Source: Notes on the practical application of… · p. 17
• The analysis window defines points contributing neighbor statistics, whereas the embedding window supplies a point’s history; thus TE can retain single-time-step resolution despite longer embedding history.Source: Notes on the practical application of… · p. 17

Repeatability of neuronal processes

• Neural contributions are generally stationary over repetitions absent learning-related drifts; split-half comparisons of means and variances can test for such drifts.The assumption concerns stability of the underlying probability distribution, not identical realizations.Source: Repeatability of neuronal processes · pp. 17–18

Relation of the ensemble method to local information dynamics

• The ensemble approach obtains local PDFs from process replications or cyclostationarity, whereas LTE localizes transfer in time given a stationary process PDF.The approaches localize different quantities and therefore can be combined.Source: Relation of the ensemble method to… · p. 18

Relation of the ensemble method to other measures of connectivity for non-stationary data

• Linear Granger causality equals transfer entropy for jointly Gaussian variables, but may miss higher-order interactions in non-Gaussian neural data; time-varying autoregressive and SETAR methods also address non-stationarity [108–111].The cited extensions fit time-dependent parameters to ensembles, while SETAR models use state-dependent parameters for transient components.Source: Relation of the ensemble method to… · pp. 18–19
• Andrzejak et al.’s repetition-based coupling measure [113] omits the target’s past and seeks a dominant coupling direction, unlike TE’s information-transfer quantification in both directions under mild conditions.TE requires the target’s past to measure information transfer; its bidirectional quantification is subject to mild conditions related to entropy production.Source: Relation of the ensemble method to… · p. 19
• The ensemble method and Strong et al.’s direct method [114] both use trials and trial-to-trial variability, but estimate different quantities: asymmetric conditional mutual information versus symmetric mutual information.Strong’s method uses intrinsic variability versus variability driven by stimulus changes.Source: Relation of the ensemble method to… · p. 19
• Strong et al.’s method [114] requires stimulus intervention, whereas intervention-free TE can quantify neural information transfer under constant stimuli; the intervention makes their method closer to causal information flow [47].Applying Strong’s method within a neural system would require precisely setting the source neuron’s or brain area’s activity.Source: Relation of the ensemble method to… · p. 19

Application of the proposed implementation to other dependency measures

• GPU nearest-neighbor searches make ensemble pooling more practical despite the high cost of generating typically about 1000 surrogate datasets, and the implementation can extend to dependency measures based on conditional mutual information.Andrzejak and colleagues [113] used ensemble pooling with delay-embedded series and nearest-neighbor statistics to estimate nonlinear dependencies.Source: Application of the proposed… · p. 19

Application to MEG source activity in a perceptual closure task

• MEG source analysis revealed time-varying information transfer, including early V1-to-OFC and V1-to-cITG connections and later exchanges between OFC/cITG and FFA.The V1-to-OFC result is consistent with a proposed OFC role in early scene segmentation [92]; cITG processes shape-from-shading information [115], and FFA is involved in face processing [116–118].Source: Application to MEG source activity in… · p. 19
• FFA appeared to be a task-network hub, receiving increasing information transfer as processing progressed, consistent with reported task-related FFA differences at about 200 ms [86].Source: Application to MEG source activity in… · p. 19
• Information-transfer patterns varied substantially across subjects, which the authors attributed to limited per-subject data; finer temporal segmentation remains a question for future work.Source: Application to MEG source activity in… · p. 19

Conclusion and further directions

• The GPU implementation of the ensemble TE method [55] reconstructed transfer in simulated and neural data over time and substantially reduced nearest-neighbor search execution times.The flexible implementation is designed to scale across experimental designs and future hardware.Source: Conclusion and further directions · p. 19
• The implementation is available for CUDA-enabled GPUs through the open-source MATLAB toolbox TRENTOOL [56].Source: Conclusion and further directions · p. 19
• The authors conclude that the ensemble method is suitable for non-stationary neural time series and may also apply to information-transfer analysis in other complex dynamic systems.Source: Conclusion and further directions · p. 20