{"paper_id":"367fd8b0-7ad6-4cf3-857a-b68f920fc972","body_text":"1 \n \n \n \n \nHeteroRC: Decoding latent information from dynamic neural responses with \ninterpretable heterogeneous reservoir computing \n \n \n \n \n \nRunhao Lu1,2 *, Sichao Liu1,3, Yanan Liu2,  \nJohn Duncan1, Richard N. Henson1,4, Alexandra Woolgar1,5 \n \n \n1 MRC Cognition and Brain Sciences Unit, University of Cambridge, Cambridge, UK  \n2 Montreal Neurological Institute, Department of Neurology and Neurosurgery, McGill University, \nMontreal, Canada \n3 Department of Production Engineering, KTH Royal Institute of Technology, Stockholm, Sweden \n4 Department of Psychiatry, University of Cambridge, Cambridge, UK \n5 Department of Psychology, University of Cambridge, Cambridge, UK \n \n* Corresponding author: Runhao Lu (runhao.lu@mrc-cbu.cam.ac.uk; rl671@cantab.ac.uk) \n  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n2 \n \nAbstract \nTime-resolved neural decoding is widely used to track information represented in neural activity, but \nconventional linear decoders primarily capture phase-locked evoked responses and often fail to recover \nrepresentations embedded in nonlinear or non–phase-locked dynamics, potentially limiting the interpretation of \nneural coding. Here, we introduce HeteroRC, a biologically inspired and interpretable decoding framework \nbased on heterogeneous reservoir computing. HeteroRC projects neural signals into a high-di mensional \nrecurrent state space with heterogeneous time constants, enabling nonlinear feature expansion and multiscale \ntemporal integration directly from raw neural time series. Simulations demonstrate that HeteroRC significantly \noutperforms linear decoders and a suite of artificial neural networks (including RNNs, LSTMs, Transformers  \nand EEGNet) on evoked responses while robustly capturing induced oscillatory power, phase synchrony, and \naperiodic modulations—dynamics that are largely latent to conventional linear methods . We further validate \nHeteroRC on two empirical EEG datasets. In a motor imagery task, it substantially improves decoding accuracy \nand exhibits superior cross-temporal generalisation, revealing dynamic representational transformations. In an \nattentional priority task, HeteroRC uncovers statistically learned spatial priority information that remains hidden \nfrom conventional methods, successfully decoding these latent states previously thought to be ‘activity-silent’. \nFurthermore, we develop a dual -level interpretability framework linking reservoir dynamics to virtual sources \nand sensor space, revealing the temporal, spectral, and spatial signatures underlying decoding performance at \nboth the individual and group levels. Together, HeteroRC offers an interpretable approach to decode information \nfrom dynamic neural responses, broadening the analytical scope of neural decoding while remaining \ncomputationally efficient and free from manual feature engineering, making it particularly suitable for small -\nsample electrophysiological studies. \n  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n3 \n \nIntroduction \nNeural decoding has become a central tool in cognitive neuroscience and brain –computer interface \nresearch, enabling inference about representational content from patterns of neural activity  [1-8] . With the \nincreasing availability of electrophysiological recordings, ranging from non-invasive magnetoencephalography \n(MEG) and electroencephalography (EEG)  to intracranial measures such as local field potentials (LFPs)  and \nmulti-unit activity , decoding methods have increasingly exploited the high temporal resolution of these \ntechniques to track neural representations over time. Time-resolved decoding and cross-temporal generalisation \nare now standard tools for characteri sing the temporal dynamics of cognitive processes such as perception, \nattention, memory, and action [6, 9-11] \nDespite their broad adoption, most time-resolved decoding studies rely on linear classifiers, typically linear \ndiscriminant analysis (LDA) or linear support vector machines (SVM) [6, 9, 12-15]. Comparative benchmarks \nhave shown that, for decoding stimulus -locked and phase-locked information expressed in evoked responses , \nlinear decoders often perform as well as or better than more complex non -linear models [16, 17] . Their \nrobustness, computational efficiency, and interpretability have therefore made linear decoding pipelines , \noperating on instantaneous signal amplitudes (e.g., voltage for EEG/LFPs or magnetic fields  for MEG), the \ndefault choice in neural time-series decoding research [6, 18]. \nHowever, this methodological standard implicitly favours a restricted class of neural signals: those that are \nphase-locked to experimental events and expressed as evoked responses [6, 15, 19]. A large body of work \nindicates that many cognitive variables are instead encoded in neural dynamics that are not phase-locked to \nstimulus onset, including induced oscillatory power [20-23], phase synchrony [24-27], and scale-free aperiodic \nactivity [28-30] . These dynamics often vary in latency and duration across trials and are therefore poorly \ncaptured by conventional  decoding approaches that operate independently at each time point on raw \ninstantaneous amplitude s using linear classifiers, which can in turn complicate the interpretation of null \ndecoding findings. In the working memory literature, for example , several studies have reported that memory \ncontent cannot be decoded from instantaneous signal traces during delay periods  [31-33]. These results have \ncontributed to the proposal that representations may be maintained in an “activity-silent” state [34], while \nremaining reactivatable by brief external perturbations (i.e., “pinging”) [31, 35, 36]. Importantly, subsequent \nwork has shown that, in some cases, information that is not decodable from raw amplitudes can be recovered \nfrom alternative signal features such as induced alpha power  [37]. T hese findings highlight that decoding \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n4 \n \noutcomes depend on the interaction between the neural signal format and the decoding model, and that failures \nof linear, instantaneous amplitude-based decoding may not uniquely determine absence of task-relevant neural \nactivity, but reflect a mismatch between representational dynamics and decoding assumptions. \nWhile non-phase-locked information can sometimes be decoded by explicitly transforming neural signals \ninto alternative features (e.g., time-frequency representations), such approaches rely on a priori assumptions \nregarding which signal dimensions carry task-relevant information. In many cognitive paradigms, however, the \noptimal representational format is unknown and may dynamically evolve across tasks, brain regions, and \nprocessing stages. Recent efforts have sought to maximi se information recovery by syste matically comparing \nand combining diverse feature sets [38, 39]. Yet, such exhaustive feature engineering remains computationally \nintensive and demonstrates that integrating multiscale features does not consistently yield additive gains  [39]. \nAlternatively, modern deep learning architectures, such as recurrent neural networks (RNNs) and Transformers, \noffer a data-driven approach to bypass manual feature extraction, utilising their high non-linear expressivity to \nlearn complex representations directly from raw time series [40]. However, the immense parameter spaces of \nthese fully trainable models inherently demand massive datasets for stable convergence. Consequently, they are \nsusceptible to overfitting and temporal smearing in the small -sample, trial-limited regimes that chara cterise \nmost cognitive electrophysiology experiments [41]. \nAddressing these limitations motivates the development of decoding frameworks that operate directly on  \nsmall-sample neural time series while remaining sensitive to non-linear, non -phase-locked, and multiscale \ntemporal dynamics. Reservoir computing (RC)  provides a biologically inspired and computationally efficient \napproach for mapping time -varying inputs into a high-dimensional dynamic state space using fixed recurrent \nconnectivity and a trained linear readout  [42-46] . Through recurrent dynamics, RC models can integrate \ninformation over time while preserving interpretability at the readout level, making them attractive for neural \ndecoding applications. However, standard RC implementations implicitly operate at a single intrinsic temporal \nscale, typically assuming homogeneous time constants across reservoir units. This assumption contrasts with \nextensive empirical and theoretical evidence that neural activity in cortex spans a hierarchy of intrinsic \ntimescales [47-49], and limits the ability of conventional RC models to represent neural dynamics that unfold \nconcurrently over fast and slow temporal scales. \nHere we introduce HeteroRC (Heterogeneous Reservoir Computing), a decoding framework that explicitly \nincorporates heterogeneous intrinsic time constants into reservoir dynamics. Motivated by empirical and \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n5 \n \ntheoretical work demonstrating a hierarchy of intrinsic timescales in cortical activity [47-50], HeteroRC samples \nreservoir units with a distribution of time constants, enabling the representation of neural dynamics across fast \nand slow temporal scales within a single recurrent state space. Multichannel neural  time series are projected \ninto this high-dimensional dynamic space, and task-relevant information is extracted using a linear readout (e.g., \nridge classification). Crucially, the linear readout enables a principled interpretation module , allowing latent \nreservoir dynamics to be linked back to s ensor-level neurophysiological signals . At the individual level, this \nframework extracts temporal, spectral, and spatial dynamics tailored to single subjects. Furthermore, by utilising \na sensor-space matching approach, these findings can be generalised to uncover shared neurophysiological \nmotifs across participants. \nThrough controlled simulations , we first demonstrate that HeteroRC significantly outperforms \nconventional linear decoders and a suite of artificial neural networks (ANNs; including RNNs, LSTMs, \nTransformers, and EEGNet) on evoked responses, while robustly capturing induced oscillatory power, phase \nsynchrony, and aperiodic modulations—neural dynamics that are typically inaccessible to conventional linear \nmethods. Applying the method to a n empirical  motor imagery neuroimaging dataset [51], we demonstrate \nimproved decoding performance relative to conventional linear decoders, particularly during internally \ngenerated imagery periods. Moreover, HeteroRC supports robust cross -temporal generalisation, capturing \nsystematic changes in representational dynamics from cue-driven to internally -maintained states. Leveraging \nthe interpretation framework, we identify the temporal, spectral, and spatial characteristics of neural signal \ncomponents that contribute most strongly to successful decoding  at the individual-subject level, and further \nreconstruct grand -average virtual sources to confirm that these distinct mechanistic signatures are robustly \nconserved across the group. Furthermore, using an attentional priority mapping dataset [36], we show that latent \ninformation not readily decodable from instantaneous amplitude-based decoding  can be recovered using \nHeteroRC, and demonstrate that decoding in perturbed and unperturbed conditions relies on partially distinct \nneural signal components, revealing different encoding regimes despite shared representational content . \nTogether, these results establish HeteroRC as a robust and interpretable framework for neural decoding, enabling \nlatent information in dynamic brain signals to be recovered and systematically interpreted beyond conventional \napproaches. \nResults \nHeteroRC \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n6 \n \nTo overcome the limitations of conventional linear decoders in capturing non-phase-locked and nonlinear \nneural dynamics, while also avoiding the requirement for large training datasets characteristic of deep learning \narchitectures, we developed HeteroRC. HeteroRC is a decoding framework that projects multichannel neural \ntime series into a high -dimensional recurrent state space governed by heterogeneous intrinsic time constants  \n(Figure 1a). These time constants are sampled from a log-normal distribution, establishing a multiscale temporal \nfilter bank that enables simultaneous integration of fast, transient responses and slower, persistent neural \ndynamics. This allows HeteroRC to operate directly on raw neural time series\n without requiring explicit feature \nengineering (e.g., pre -computed spectral power or phase), while retaining the computational efficiency of a \nfixed reservoir coupled with a trained linear readout. \nTo address the interpretability challenges associated with recurrent models, we developed an interpretation \nframework that links  reservoir dynamics to physiological signal generators (Figure  1b). First, we apply a \ncovariance-based activation pattern analysis [18] to project the learned readout weights back into reservoir space \nand identify units that contribute most strongly to decoding.  By clustering the temporal dynamics of these \ninformative units, we extract latent virtual source signals that capture task -relevant dynamics within the \nreservoir. These virtual sources can be analysed directly in the time and frequency domains and are further \nprojected back to sensor space to reveal their corresponding spatial topographies . Together, this framework \ndisentangles the temporal, spectral, and spatial signatures underlying decoding performance, allowing phase -\nlocked evoked responses to be distinguished from induced oscillatory and aperiodic components, and directly \nlinking decoding outcomes to interpretable neural signal mechanisms. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n7 \n \n \nFigure 1. HeteroRC decoding and interpretation framework. \n(a) Decoding framework. Multichannel neural signals (trials × channels × time) are provided as inputs to a \nrecurrent reservoir. Inputs are linearly projected to reservoir units through a fixed, randomly initialis ed input \nweight matrix ( 𝑊𝑊in ), with full connectivity between input channels and reservoir units. Reservoir units are \nconnected via a fixed recurrent weight matrix ( 𝑊𝑊res ) with sparse random connectivity (10% non-zero \nconnections), scaled to ensure stable echo-state dynamics. Each reservoir unit is endowed with an intrinsic time \nconstant (𝜏𝜏), sampled from a log-normal distribution, resulting in heterogeneous temporal integration properties \nspanning fast and slow timescales. Neural inputs are thus transformed into high -dimensional reservoir state \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n8 \n \ntrajectories that capture multiscale temporal dynamics. At each time point, reservoir states are decoded using a \ntrained linear readout ( 𝑊𝑊out ; ridge classification), yielding time -resolved decoding performance for task \nconditions or stimulus classes.  (b) Interpretation framework. Trained readout weights are first projected back \ninto reservoir space using a covariance-based activation pattern analysis (Haufe’s transform), identifying \nreservoir units that contribute most strongly to decoding. From here, the framework branches into two analytical \ntracks. For individual-level interpretation, t he dynamics of informative units are directly clustered to extract \nindividual-specific latent virtual source signals. These virtual sources are analysed  in the time and frequency \ndomains (e.g., evoked responses, time -frequency representations, fitting oscillations & one over  (FOOOF) \nanalysis) and further projected back to sensor space to reveal their corresponding spatial topographies, enabling \njoint temporal, spectral, and spatial interpretation of decoding performance.  For group-level interpretation, \nselected units from all participants are first projected to a common sensor space. These spatial topographies are \npooled and globally clustered, and the resulting cluster assignments are mapped back to individual reservoirs to \nreconstruct group-level virtual source signals.\n These sources are subsequently subjected to the same temporal, \nspectral, and spatial analyses mentioned above. \nHeteroRC decodes diverse classes of neural dynamics under controlled simulations \nTo systematically evaluate the sensitivity of HeteroRC to distinct classes of neural dynamics, we generated \nsynthetic datasets simulating five canonical  signal types commonly observed in neural time series recordings: \n1) phase-locked evoked responses, 2) induced oscillatory power, 3) inter -site phase clustering (ISPC), and 4) \nslope and 5) offset of aperiodic spectral modulations . Simulated signals were embedded in realistic 1/ f \nbackground activity with white noise over 0.8-s epochs, with task-relevant modulations confined to a 0.2–0.6-\ns time window and subject to trial -by-trial temporal jitter (see Methods). For simulations involving oscillatory \npower and phase synchrony, task -relevant modulations were centred at 10 Hz, reflecting the canonical alpha-\nband activity ; however, identical decoding patterns were obtained when modulations were placed at other \nfrequencies (Figure S1) . Crucially, the strength of each signal type was manipulated independently while \nholding the remaining components constant, ensuring that decoding performance reflected sensitivity to the \ntargeted dynamical feature.  We employed an event-related design with 5 -fold cross-validation, where feature \nscaling and model training were strictly separated within each fold to prevent data leakage.  \nWe first evaluated decoding performance in simulations where discriminative information was carried by  \nphase-locked evoked responses , a regime typically considered the optimal use case for linear , instantaneous \namplitude-based classifiers. In this setting, HeteroRC outperformed both LDA and SVM , achieving higher \ndecoding accuracy within the task-relevant window while preserving accurate temporal locali sation of the \nevoked response peak (Figure 2a). \nA pronounced divergence in performance between HeteroRC, LDA and SVM emerged when \ndiscriminative information was embedded in non-phase-locked dynamics, including induced oscillatory power \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n9 \n \nwith randomised phase, ISPC, and aperiodic slope and intercept modulations. Under these conditions, both LDA \nand SVM applied to raw time series failed to recover task -relevant information, yielding performance near \nchance levels (Figure 2b–d). In contrast, HeteroRC robustly decoded these latent dynamics, despite the absence \nof consistent phase-locked amplitude differences. It moreover accurately recovered the time-window in which \nthe latent dynamics had been applied. \nTo explicitly contrast HeteroRC with conventional feature engineering pipelines, we additionally evaluated \nthe performance of LDA when it was trained and tested on data that was manually  transformed into time -\nfrequency space – specifically, on 8–12 Hz alpha-band power extracted via Morlet wavelets and a Hilbert filter \napproach (Figure S2). We again tested this approach for detecting underlying effects embedded in evoked \nresponses, induced oscillatory power with randomi sed phase, ISP C, and aperiodic  slope and intercept \nmodulations. As expected, the feature engineering approach successfully decoded information when the \nunderlying neural dynamics  matched the features of the data that the manual engineering approach targeted . \nThus, LDA based on 8-12Hz alpha-band power was able to discriminate the classes when the underlying signal \nvaried in 10 Hz induced power, and when the underlying modulation was a change in the broadband aperiodic \nintercept, since this manipulation also  altered absolute alpha-band power . It could also weakly and partially \ndetect the underlying modulation of aperiodic slope, for the same reason. However, the engineering pipeline  \ncompletely failed to recover phase-locked evoked responses, ISPC, and induced power outside the pre-specified \nfilter range (e.g., 20 Hz). Furthermore, even when decoding was successful (e.g., for 10 Hz power), both time-\nfrequency methods, and particularly the wavelet approach, exhibited pronounced temporal smearing extending \nbeyond the ground-truth 0.2–0.6-s window. These results demonstrate that explicit feature extraction is not only \nbottlenecked by manual parameterisation but also susceptible to integration-induced temporal blurring. \nGiven this observation, and because HeteroRC  also inherently relies on its own recurrent temporal \nintegration, an important concern is whether its robust decoding of latent dynamics could be similarly \nconfounded by temporal smearing or systematic latency shifts. To address this, we compared different reservoir \nprocessing strategies in controlled simulations, including unidirectional processing and alternative bidirectional \nfusion schemes (see Methods for details).\n As shown in Figure S3, in simulations dominated by evoked responses \nand aperiodic activity, unidirectional reservoirs exhibited systematic temporal lags and broadened decoding \npeaks, consistent with integration-induced smearing. Bidirectional processing substantially reduced these \neffects by cancelling  direction-specific temporal biases. Among bidirectional strategies, multiplicative fusion \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n10 \n \n(the method used in Figure 2) provided sharper temporal localisation than averaging (as in Figure S3), while \npreserving robust decoding performance. These results indicate that, unlike conventional filtering approaches, \nthe sustained decoding observed with HeteroRC more faithfully reflects true underlying neural dynamics rather \nthan algorithmic lag, thereby motivating the use of bidirectional multiplicative fusion throughout the main \nanalyses. \nTogether, these simulations demonstrate that HeteroRC not only exceeds the performance of conventional \nlinear decoders in their preferred regime of phase-locked evoked responses, but it also generalizes to reliably \nrecover information encoded in non-phase-locked and aperiodic neural dynamics. This broad sensitivity enables \naccurate, time-resolved decoding across a wide range of biologically meaningful signal formats commonly \npresent in neural time series recordings. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n11 \n \n \nFigure 2. HeteroRC decodes diverse classes of neural dynamics in controlled simulations \nSimulated datasets were generated to isolate and test decoding sensitivity to five distinct classes of neural \ndynamics: phase-locked evoked responses (a), induced oscillatory power with randomi sed phase (b), inter-site \nphase clustering (ISPC) (c), and aperiodic spectral modulations affecting the 1/ f slope (d) and offset (e). In all \nconditions, task-relevant modulations were confined to a ~0.2 –0.6-s time window and subject to trial -by-trial \ntemporal jitter. For simulations involving oscillatory power and ISPC, task -relevant modulations were centred \nat 10 Hz (for results at different frequencies see Figure S 1). For each simulation, the left panels illustrate the \nunderlying signal differences between conditions, including evoked responses (time domain), power spectral \ndensity (PSD, frequency domain), and ISPC  (phase synchrony), while polar plots depict phase distributions. \nRight panels show time -resolved decoding accuracy for HeteroRC (red) compared with conventional linear \ndecoders (LDA, blue; SVM, purple). Horizontal bars indicate time points with decoding performance \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n12 \n \nsignificantly above chance (p < 0.05, corrected by cluster-based permutation test).  \nHeteroRC outperforms artificial neural networks in decoding evoked responses \nWhile ANNs offer high nonlinear expressivity, their effectiveness inherently scales with the availability of \nlarge training datasets. To contextualise HeteroRC against this class of machine learning approaches, we \ncompared its performance against four canonical ANNs: RNNs [52] and Long Short -Term Memory (LSTM) \nnetworks [53], which rely on recurrent hidden states to integrate temporal sequences; Transformers [54], which \nuse self-attention mechanisms capture global temporal dependencies; and EEGNet [55], a convolutional neural \nnetwork (CNN)-based architecture specifically optimised for the spatial and temporal constraints of EEG signals. \nWe focused this comparative benchmark specifically on phase- locked evoked responses , as these highly \nconsistent signals represent the most straightforward scenario for ANNs to achieve stable performance. Due to \nits convolutional architecture, EEGNet was evaluated exclusively using a windowed formulation to generate a \ntime-resolved decoding profile comparable to HeteroRC. In contrast, the other models were assessed in both \nstandard sequence-to-sequence and windowed regimes. \nFor the general-purpose ANNs (RNN, LSTM, and Transformer), we first evaluated them in a standard \nsequence-to-sequence regime. In this approach, the models process the full temporal epoch at once and are \ntrained to output a continuous trajectory of class pr edictions, generating a discrete decision at every individual \ntime step . For the RNN and LSTM, this prediction relies on the accumulated signal history, whereas the \nTransformer globally integrates both past and future context across the entire epoch.  In this setting, HeteroRC \nsubstantially outperformed the se models (Figure 3 a). Notably, RNN and LSTM  exhibited a pronounced \ntemporal lag and smearing, with their decoding accuracy peaking later but persisting longer than the ground-\ntruth neural modulation. This reflects the integration time required to accumulate evidence from past inputs and \nthe persistence of information within their hidden states even after the underlying signal had ceased. The \nstandard Transformer poorly localized neural dynamics in time, as its unmasked global self-attention mechanism \nintegrates both past and future context across the entire epoch. \nTo mitigate these issues and provide a more rigorous baseline, we next trained windowed variants of the \nANNs alongside EEGNet. By applying a 100ms sliding window across the epoch, we strictly bounded the \ntemporal integration horizon, forcing the models to rely on local signal features rather than accumulated history \nor future context. Although this windowing strategy improved the temporal localization of the decoding results,\n \npartially correcting  the temporal lag and smearing , HeteroRC continued to demonstrate superior decoding \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n13 \n \naccuracy and temporal precision (Figure 3b). \nThese results indicate that HeteroRC, utilising a fixed reservoir with heterogeneous intrinsic time constants, \ncaptures multiscale neural dynamics more efficiently than fully trainable ANNs when data is limited . By \nintegrating these dynamics with bidirectional multiplicative fusion to actively cancel temporal lags, HeteroRC \nachieves a highly favourable trade-off between representational richness, data efficiency, and temporal precision. \n \nFigure 3. HeteroRC outperforms ANNs in decoding evoked responses \nTime-resolved decoding accuracy for simulated phase -locked evoked responses, comparing HeteroRC (red) \nagainst standard and windowed ANNs. Task-relevant modulations were confined to a ~0.2–0.6-s time window. \n(a) Comparison with standard sequence-to-sequence models: Recurrent Neural Network (Std RNN, blue), Long \nShort-Term Memory network (Std LSTM, purple), and Transformer (Std Transformer, green). (b) Comparison \nwith windowed (sequence-to-one) models utilising a 100-ms sliding window, including Windowed RNN (blue), \nLSTM (purple), Transformer (green), and the EEGNet (yellow).  EEGNet was evaluated exclusively using the  \nwindowed formulation to respect its convolutional architecture. Shaded regions denote the standard error of the \nmean across simulated subjects. Horizontal dashed lines indicate chance-level performance (0.5). Solid \nhorizontal bars at the bottom indicate time windows with robust decoding performance for the correspondi ng \nmodels.  \nHeteroRC decodes sustained internally generated motor imagery representations \nTo validate HeteroRC on empirical data, we first applied the framework to a motor imagery dataset ( N = \n9; BCI Competition IV 2a  [51]), which requires decoding of four classes of imagined movements (left hand, \nright hand, feet, and tongue) from a 22-channel EEG (Figure 4a). Motor imagery relies on internally generated \nneural activity with variable onset and duration across trials, providing a stringent test for decoding non-phase-\nlocked dynamics. In addition, training and evaluation data were collected in separate recording sessions on \ndifferent days for each participant, requiring models to generali se across sessions despite potential non -\nstationarities in signal quality or electrode impedance.\n Given that our controlled simulations demonstrated near-\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n14 \n \nidentical performance between LDA and linear SVM, we restricted our subsequent comparisons to LDA for \ncomputational efficiency. \nWe first evaluated time-resolved decoding accuracy, comparing HeteroRC against LDA (Figure 4b). In the \nexperimental paradigm, each trial consisted of a fixation period, followed by a visual cue indicating the \nmovement to be imagined , and a sustained motor imagery period in the absence of external sensory input. \nFollowing cue onset, both HeteroRC and LDA showed a rapid increase in decoding accuracy with comparable \nlatencies. However, LDA performance declined to near-baseline levels before cue offset, whereas HeteroRC \nmaintained robust and significant decoding accuracy. Crucially, during the sustained imagery period (from 1.25 \ns\n onward), when no external stimulus was present and decoding depended entirely on internally maintained \nrepresentations, HeteroRC significantly outperformed the linear decoder.  HeteroRC sustained above-chance \ndecoding for more than 1.5 s, whereas LDA exhibited reliable decoding for only approximately 250 ms. In \naddition, inspection of individual-subject decoding dynamics (Figure S4) suggested that HeteroRC numerically \noutperformed LDA in decoding accuracy for all participants. To quantify these individual-subject level effects, \nwe compared peak decoding accuracy between HeteroRC and LDA, separately for the cue period (0 –1.25 s) \nand the sustained imagery period (1.25 –3 s). HeteroRC showed significantly higher peak decoding accuracy \nthan LDA in both cue (t(8)= 5.32, p < 0.001) and imagery periods (t(8) = 3.91, p = 0.004) (Figure 4b, right). \nThis indicates that the sustained decoding advantage of HeteroRC is robust across participants and not driven \nby a small subset of subjects. \nTo further characterise the temporal structure of the underlying neural codes, we computed cross-temporal \ngeneralisation matrices ( Figure. 3c). The linear decoder exhibited a predominantly diagonal generali sation \npattern confined to the cue period and the early imagery phase, indicating reliance on transient , time-specific \nand likely phase -locked features. In contrast, HeteroRC revealed a dynamic evolution in representational \nstructure. During the early cue phase (before ~600 ms), generalisation was largely diagonal, consistent with \ndynamically evolving sensory representations. This pattern progressively transitioned into a broad, square-like \ngeneralisation structure during the late cue phase (after ~600 ms) and persisted throughout the sustained imagery \nperiod, indicating the emergence of a temporally stable representational regime.  \nTo identify the neural mechanisms supporting this sustained decoding during imagery, we applied the \ninterpretation framework to an individual  participant (Subject 1) at the peak decoding time point  within the \nimagery period (2.3 s). Using a covariance-based activation pattern analysis [18], we identified the top 25 most \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n15 \n \ninformative reservoir units and grouped them into three functional clusters via hierarchical clustering, yielding \nthree virtual source signals that captured distinct latent dynamics contributing to decoding (see Methods for \ndetails). Inspection of these virtual sources revealed that HeteroRC leveraged multiple signal characteristics to \ndistinguish imagined motor commands. Some components (e.g., Cluster 1) exhibited distinct condition-specific \ntemporal modulations, whereas some reflected broader spectral dynamics (Cluster 2). Notably, Cluster 3 \ncaptured differences in oscillatory power within the alpha/mu range (8-13 Hz) as well as class-dependent \nvariations in aperiodic spectral slope, which differentiated  hand imagery from foot and tongue imagery .\n \nProjection of these sources back to sensor space further revealed spatially distinct cortical topographies \nassociated with each d ynamical component. Together, these results indicate that the reservoir integrates \ninformation across temporal, spectral, and spatial domains to form a robust representation of internally \ngenerated motor states. \nIn summary, HeteroRC robustly decodes internally generated motor imagery representations that are \ndifficult to recover using conventional linear models. Beyond improved decoding accuracy, cross -temporal \ngeneralisation analyses demonstrate that HeteroRC tracks the transformation of neural representations from \ntransient cue-driven activity to stable, internally maintained imagery states. The accompanying interpretation \nframework further reveals the distinct temporal, spectral, and spatial signal components supporting thi s \nperformance, establishing HeteroRC as a powerful and interpretable tool for characterising cognitive processes \nthat unfold in the absence of external sensory input. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n16 \n \n \nFigure 4. HeteroRC decodes internally generated motor imagery and captures representational dynamics \n(a) Experimental design of the motor imagery task [51]. Each trial consisted of a fixation period, followed by a \nvisual cue (an arrow pointing left, right, down, or up) indicating the movement to be imagined (left hand (LH), \nright hand (RH), feet, or tongue, respectively), and a sustained motor imagery period in the absence of external \nsensory input. (b) Left: Time-resolved decoding accuracy over time for HeteroRC (red) and linear discriminant \nanalysis (LDA; blue). Shaded regions denote the standard error of the mean across participants. Horizontal bars \nindicate time points with decoding accuracy significantly above chance ( p < 0.05, corrected  by cluster-based \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n17 \n \npermutation test) for HeteroRC (red), LDA (blue) and the difference between them (black) . Right: Peak \ndecoding accuracy for each participant during the cue period and the sustained imagery period. Each dot \nrepresents one participant, with lines connecting LDA and HeteroRC within subjects. HeteroRC shows \nsignificantly higher peak decoding accuracy than LDA in both periods (\n**p < 0.01; ***p < 0.001). (c) Cross-\ntemporal generalisation matrices. Decoding accuracy as a function of training time (y-axis) and testing time (x-\naxis) for LDA (left), HeteroRC (middle), and their difference (HeteroRC minus LDA, right). Black contours \nindicate regions significantly above chance  (p  < 0.05, corrected  by cluster -based permutation test ). (d) \nInterpretation of reservoir dynamics during sustained imagery for a n individual  participant. Top: Haufe-\ntransformed activation patterns plotted  as a function of reservoir unit intrinsic timescale, with most 25 \ninformative units grouped into three clusters using hierarchical clustering.  Bottom: Analyses of virtual source \nsignals derived for each cluster, showing trial-averaged temporal profiles, time– frequency representations, \npower spectral density estimates of oscillatory and aperiodic components (derived using FOOOF), and \ncorresponding sensor-space projections. \n \nHeteroRC uncovered latent attentional priority representations invisible to linear decoding  \nFinally, we evaluated whether HeteroRC could recover task-relevant information from neural states that \nare typically considered inaccessible to standard raw amplitude-based decoding approaches. We addressed this \nquestion using an attentional priority mapping dataset  (N = 23, Figure 5a) [36] in which participants acquired \nspatial priority maps through statistical learning. Briefly, participants learnt implicitly that one of eight locations \nin a visual display  was more likely to contain a target (a singleton shape).  Previous analyses using linear \nclassifiers reported that neural representations of these priority locations was undetectable during the inter-trial \ninterval and could only be detecting following a brief visual impulse (“ping”), leading to the proposal that the \nrepresentations were maintained in an activity-silent state. \nWe compared the decoding of spatial priority between HeteroRC and LDA  in both Ping and No-Ping \nconditions (Figure 5b). Consistent with the original publication [36], LDA successfully decoded spatial priority \nfollowing the visual impulse but failed to achieve above -chance performance in the absence of the impulse, \nsuggesting the absence of decodable information during No-Ping trials .\n In contrast, using HeteroRC, spatial \npriority information was robustly decoded in both Ping and No-Ping conditions. There was significant decoding \naccuracy observed throughout the epoch in both cases, including before the ping onset (and equivalent no-ping \ntimepoint). Nonetheless, the ping still elicited a transient numerical increase in decoding accuracy. \nAt the individual-subject level, again, HeteroRC tended to show higher decoding accuracy than LDA for \nnearly all participants in both Ping and No-Ping conditions, as illustrated by per-subject time-resolved decoding \nresults (Figure S5). To quantify these effects across participants, we compared peak decoding accuracy between \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n18 \n \nHeteroRC and LDA separately for Ping (t(22)  = 8.92, p < 0.001) and No-Ping (t(22) = 6.28, p < 0.001) trials. \nHeteroRC showed significantly higher peak decoding accuracy than LDA in Ping and No-Ping conditions \n(Figure 5c). \nTo identify the neural features supporting this latent information, we applied the interpretation framework \nto an individual participant (Subject 24) at the peak decoding time point for Ping (0.34 s) and No-Ping (0.48 s) \ntrials. In Ping trials (Figure 5d), virtual source analysis revealed prominent evoked responses that discriminated \npriority locations (both clusters), consistent with phase-locked reactivation of the priority map by the external \nimpulse. This involved spatially distributed generators across frontal and posterior regions. In contrast, No-Ping \ntrials (Figure 5e , both clusters ) showed no discernible evoked components. Instead, decoding was driven by \ninduced neural dynamics: virtual sources exhibited condition-specific oscillatory activity in the alpha (~10 Hz) \nand beta (~20 Hz) bands, accompanied by shifts in aperiodic spectral slope and offset. Projection of these \nsources to sensor space localised the dominant contributions primarily to posterior cortical regions. \nThese results suggest that neural information traditionally considered inaccessible to instantaneous \namplitude-based linear decoding (e.g., attentional priority information) may in fact be continuously maintained \nin latent, non -phase-locked dynamics, which can be readily recovered by HeteroRC.  On the other hand, t he \naccompanying interpretation framework further reveals that decoding in Ping and No-Ping conditions is \nsupported by distinct neural signal components. This observation adds important new insight to the debate \nsurrounding activity-silent mechanisms, suggesting both that previously “hidden” information may in fact have \nbeen maintained in non-phase locked activity, and also that additional evoked activity, more closely aligned \nwith activity-silent mechanisms, can be recovered through pinging. These observations showcase the additional \ninsight possible through the proposed interpretability approach using HeteroRC. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n19 \n \n \nFigure 5. HeteroRC recovers latent task information in an attentional priority mapping task. \n(a) Experimental design of the attentional priority mapping task  [36]. Participants had to report the orientation \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n20 \n \nof a line surrounded by a singleton shape. Participants implicitly learned spatial regularities in target locations \nacross blocks, forming a learned attentional priority map. In Ping trials, a brief visual impulse (200 ms) was \npresented during the inter -trial interval. In No -Ping trials, no visual stimulus was presented, but a matched \ntrigger event was recorded at an equivalent time point .\n The subsequent visual search display required \nparticipants to report the target feature and was not used for decoding. (b) Time-resolved decoding accuracy of \nspatial priority for linear discriminant analysis (LDA; top) and HeteroRC (bottom) in Ping (red) and No -Ping \n(grey dashed) conditions, aligned to ping (or no-ping interval)  onset. Shaded regions denote the standard error \nof the mean across participants. Horizontal bars indicate time points with decoding accuracy significantly above \nchance (p < 0.05, corrected by cluster-based permutation test). The LDA results replicate the results reported in \nthe original paper, in which spatial priority can only be decoded in the ping (red) condition. However, HeteroRC \nshows sustained decoding of attentional priority before and throughout this time window, for both ping (red) \nand no-ping (grey). (c) Peak decoding accuracy for each participant in Ping (top) and No-Ping (bottom) \nconditions. Each dot represents one participant, with lines connecting LDA and HeteroRC within subjects. \nHeteroRC shows significantly higher peak decoding accuracy than LDA in both conditions (\n***p < 0.001). (d) \nInterpretation of reservoir dynamics in Ping trials  for a n individual  participant. Top: Haufe -transformed \nactivation patterns plotted as a function of reservoir unit intrinsic timescale, with most 25 informative units \ngrouped into two clusters using hierarchical clustering. Bottom: Analyses of virtual source signals derived for \neach cluster, showing trial-averaged temporal profiles, time–frequency representations, power spectral density \nestimates of oscillatory and aperiodic components (derived using FOOOF), and corresponding sensor-space \nprojections. These clusters primarily capture distinct phase-locked evoked responses . (e) Interpretation of \nreservoir dynamics in No -Ping trials  (same analyses as in d). In contrast, these clusters capture non-phase -\nlocked induced alpha/beta dynamics and aperiodic shifts. \n \nGroup-level interpretation of latent dynamics via sensor-space matching \nWhile the single -subject interpretation above  provides highly resolved spatiotemporal signatures, a \nfundamental challenge in applying reservoir computing to neural data is generalising these mechanistic motifs \nacross groups of participants. Because the randomly initialised latent reservoir spaces are idiosyncratic to each \nparticipant, direct cross -subject comparison of internal reservoir units is challenging . To overcome this, we \ndeveloped a sensor -space matching approach. We extracted the top 10 most informative units from each \nparticipant and projected their activity back to the standardi sed EEG sensor level, a physical space that is \nanatomically comparable across individuals. We then performed a global clustering analysis on the spatial \ntopographies of all extracted units, grouping them into distinct spatial clusters  representing shared \nneurophysiological motifs. \nBy mapping these spatially matched units back to their respective reservoirs, we reconstructed grand -\naverage virtual source signals across participants (Figure 6). Applying this framework first to the motor imagery \ndataset, the group-level interpretation isolated representational motifs across the cohort that complemented the \nhighly specific dynamics observed at the individual level. Specifically, the global spatial clustering disentangled \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n21 \n \nthree distinct topographical components: two laterali sed clusters peaking over the left (Cluster 1) and right \n(Cluster 3) central electrodes, and a centrally distributed cluster over the midline (Cluster 2). These distinct \nspatial topographies cleanly indicated  lateralised hand imagery from centrally represented foot and tongue \nimagery. Furthermore, analysis of the reconstructed virtual sources revealed that these spatial clusters captured \nfunctionally distinct dynamical profiles. Consistent with the single-subject observations, some components (e.g., \nCluster 2) exhibited more pronounced condition-specific temporal modulations, whereas the laterali sed \ncomponents (Clusters 1 and 3) were characterised by strong, oscillatory power modulation in the alpha/mu and \nbeta bands.  \nSimilarly, in the attentional priority dataset, the group -level virtual sources corroborated the mechanistic \ndivergence between maintenance states  seen in the example individual subject level above . Decoding in the \nPing condition was predominantly driven by strong, phase -locked evoked responses locali sed to posterior \nsensors (Figure 6b ). In contrast, decoding in the No-Ping condition reflected  induced oscillatory power \nmodulations (Figure 6c).  \nCollectively, these findings not only reinforce that successful decoding of neural representations relies on \nthe joint contribution of diverse temporal, spectral, and spatial neural signatures , but also demonstrate that \nHeteroRC can reliably extract and disentangle heterogeneous neural codes across individuals without requiring \ndirect hyperalignment of their idiosyncratic latent state spaces. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n22 \n \n \nFigure 6. Group-level interpretation of latent dynamics via sensor-space matching. \n(a) Experimental Group-level interpretation of reservoir dynamics in the motor imagery dataset across all \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n23 \n \nparticipants. Top: Haufe-transformed activation patterns of the top 10 most informative units per participant \nplotted as a function of reservoir unit intrinsic timescale, grouped into three global clusters based on their \nspatial topographies. Bottom: Analyses of grand-average virtual source signals reconstructed for each spatial \ncluster, showing trial -averaged temporal profiles, time– frequency representations, power spectral density \nestimates of oscillatory and aperiodic components (derived using FOOOF), and corresponding grand-average \nsensor-space projections.\n (b) Group -level interpretation of reservoir dynamics in Ping trials from the \nattentional priority dataset (same analyses as in a, grouped into two global clusters). (c) Group -level \ninterpretation of reservoir dynamics in No-Ping trials. \nDiscussion \nNeural decoding is widely used to infer when and how information is represented in the brain, particularly \nin time -resolved analyses of EEG, MEG, and LFPs signals [2, 6, 9] . Conventional decoding pipelines in \nneuroscience predominantly rely on linear classifiers applied to sensor -level raw amplitude time series, which \nare effective for capturing stimulus-locked activity [16, 17] but are less sensitive to information encoded in non-\nphase-locked or nonlinear neural dynamics  [21, 37]. As a consequence, a failure to decode information using \nstandard pipelines risks being misinterpreted as an absence of active neural coding. Our findings caution against \nthis, demonstrating that decodability is not solely a property of the neural representation itself, but depends \ncritically on the interaction between neural dynamics and the assumptions of the decoding model. Across \ncontrolled simulations and two empirical datasets, HeteroRC consistently recovered task -relevant information \nthat was inaccessible to conventional linear amplitude -based decoders, including information encoded in \ninduced oscillatory activity, phase synchronisation, aperiodic dynamics, and internally generated neural states.  \nCrucially, the accompanying individual- and group-level interpretation framework enabled these decoding \nresults to be linked back to identifiable temporal, spectral, and spatial neural signatures, revealing distinct \ndynamical regimes supporting decoding under different task conditions. \nA key feature of the present decoding framework is its ability to operate directly on raw neural time-series \ndata without explicit feature engineering, while remaining sensitive to information expressed in both phase-\nlocked and non-phase-locked and nonlinear neural dynamics. Rather than predefining representational domains \n(e.g., power, phase, or connectivity), HeteroRC embeds neural time series signals into a recurrent state space \nwith fixed random connectivity and heterogeneous intrinsic time constants. The use of heterogeneous intrinsic \ntime constants is motivated by a growing body of empirical and theoretical work demonstrating that neural \ndynamics in the cortex unfold across multiple, partially overlapping timescales. Single-neuron recordings have \nrevealed substantial variability in intrinsic and effective time constants, supporting the idea that cortical circuits \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n24 \n \nact as a reservoir of temporal integration windows capable of maintaining information over different durations  \n[50]. At the population and systems level, recent work has shown that intrinsic timescales are organi sed \nhierarchically across cortical areas, with progressively longer integration timescales in higher -order regions, \nfacilitating reliable signal propagation and accumulation of task -relevant information [47, 49] . Theoretical \nanalyses further suggest that heterogeneity in neuronal t ime constants is not a nuisance but a computational \nresource that shapes population dynamics and expands the repertoire of representational transformations \navailable to neural circuits [48]. By incorporating a distribution of intrinsic time constants, HeteroRC introduces \na biologically motivated inductive bias that reflects multiscale cortical dynamics, allowing information \ndistributed over fast and slow temporal regimes to be accessed by linear decoding without assuming \ninstantaneous signal expression. \nA known limitation of recurrent decoding approaches is that temporal integration of past inputs can \nintroduce systematic temporal smearing, reducing the precision with which the timing of informative neural \nevents can be recovered [42]. To mitigate this effect\n in offline analyses, HeteroRC incorporates bidirectional \ntemporal processing, in which neural time series are processed both forward and backward in time, and the \nresulting reservoir states are combined to cancel direction -specific temporal lags. We explicitly note that this \nbidirectional, multiplicative fusion is a non -causal technical strategy rather than a biologically plausible \nmechanism. Conceptually analogous to zero-phase filtering in standard signal processing  [56], it is employed \nhere to enhance temporal resolution during offline decoding while retaining the expressive benefits of recurrent \nintegration. For applications requiring strict biological plausibility or real-time, online decoding (e.g., brain –\ncomputer interfaces), the framework readily defaults to a causal, unidirectional approach. Ultimately, by training \nonly a regulari sed linear readout, the framework balances expressive power with data efficiency and \ninterpretability [18], making it well suited to the small-sample regimes and inferential goals typical of cognitive \nneuroscience experiments. \nThe decoding results point to a common principle: task-relevant information in neural time-series is not \nuniformly expressed as transient, phase-locked responses, but often emerges through sustained, trial -variable \nneural dynamics. Linear classifiers applied to sensor-level amplitude traces are well suited to decoding evoked \nresponses, as evidenced by both simulations and empirical data, where linear decoding captured early, stimulus-\nlocked information with high temporal precision [14, 16, 17]. However, beyond this initial evoked regime, \ndecoding performance with linear models rapidly declined, suggesting limited sensitivity to neural information \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n25 \n \ndistributed over time or expressed in non-phase-locked dynamics. In contrast, HeteroRC consistently recovered \ntask-relevant information across a broader range of neural regimes. In the motor imagery dataset, HeteroRC \nmaintained robust decoding throughout extended imagery periods in the absence of external stimulation, \nconsistent with the view that internally generated motor representations are supported by sustained and induced \nneural dynamics rather than brief cue-locked responses [7, 19, 57]. Similarly, in the attentional priority mapping \ndataset, HeteroRC decoded spatial priority information both following a visual impulse and during no-ping \nintervals, indicating that learned priority representations remained accessible even when they were n ot \nexpressed as overt phase-locked responses. These convergent findings across tasks suggest that many neural \nrepresentations previously considered difficult to decode may persist as continuous, dynamically evolving states \nthat are poorly matched to the assumptions of instantaneous amplitude -based linear decoding.\n Together, the \nsimulation and empirical results indicate that the apparent boundaries of decodability in neural time -series \nsignals are shaped not only by the presence or absence of neural information, but by how that information is \ndynamically expressed and sampled by the decoding model. By integrating neural activity over time while \npreserving temporal precision, HeteroRC provides access to representational formats that extend beyond evoked \nresponses and are central to internally generated and learned cognitive states. \nBeyond decoding performance, an important contribution of the present work lies in the interpretability \nframework that links successful decoding to identifiable neural signal components. In the context of \nEEG/MEG/LFPs decoding, improved accuracy alone provides limited insight into how information is \nrepresented in neural activity . By projecting decoding weights back into reservoir space and further mapping \ninformative reservoir dynamics to latent virtual sources and sensor -level patterns, the present framew ork \nconstrains decoding results to be interpretable in terms of temporal, spectral, and spatial neural signatures.\n \nApplying this framework revealed that decoding success in different task contexts relied on distinct classes of \nneural dynamics. In the motor imagery dataset, informative components reflected a combination of sustained \ntemporal dynamics, oscillatory power modulations in sensorimotor rhythms, and aperiodic spectral changes. In \nthe attentional priority mapping dataset, interpretation distinguishe d between evoked, phase -locked responses \nfollowing an external impulse and induced, non-phase-locked dynamics supporting decoding in the absence of \nperturbation. The fact that Ping and No -Ping decoding relied on partially distinct neural signal components \nprovides convergent evidence that HeteroRC accesses multiple representational regimes rather than exploiting \na single dominant feature. By extending this interpretation framework to the group level via sensor-space spatial \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n26 \n \nmatching, we further demonstrated that these distinct mechanistic motifs are conserved across the cohort , \nreflecting a generalised neurophysiological phenomenon rather than an idiosyncratic, single -subject artifact. \nMore broadly, these results illustrate how interpretability can serve as a critical safeguard for decoding -based \ninference in neuroscience. By revealing which neural dynamics support decoding under different conditions, \nboth at the individual and group levels , the interpretation framework enables researchers to assess whether \ndecoding results are consistent with known physiological mechanisms and task demands. In this way, \ninterpretability is not merely a descriptive add -on, but an essential component for drawing meaningful \nconclusions about neural representations from time-resolved decoding analyses. \nWhile recent advances in neural decoding have increasingly explored large-scale nonlinear models, \nincluding CNNs, RNNs,  and transformer -based architectures, these approaches are typically optimi sed for \nepoch-wise classification and often obscure fine-grained temporal correspondence through hierarchical \nconvolution and pooling operations, complicating analyses of representational dynamics and cross -temporal \ngeneralisation\n [40, 55, 58] . Moreover, their effectiveness commonly depends on large training datasets and \nextensive hyperparameter optimisation, which are frequently incompatible with the small-sample regimes, inter-\nsubject variability, and interpretive goals characteristic of neuroscience experiments [40, 41]. Our comparative \nsimulations directly substantiate these concerns.  Specifically, we showed that standard sequence-to -sequence \nANN architectures are highly susceptible to pronounced temporal lags. Because recurrent networks (such as \nLSTMs) continuously accumulate all preceding information, and unmasked Transformers globally integrate \nboth past and future context, their resulting decoding time courses inherently misalign with the ground-truth \ntemporal windows of the underlying neural signals. Although restricting the temporal integration horizon via a \nsliding-window approach successfully optimises the temporal precision of these trainable models, their overall \ndecoding accuracy and temporal sharpness still underperform relative to HeteroRC.\n In comparison to these \nmodels, HeteroRC occupies a complementary position in the modelling landscape. By relying on fixed random \nrecurrent connectivity and training only a linear readout, HeteroRC avoids large -scale parameter optimisation \nwhile still enabling rich temporal representations through recurrent dynamics. This design yields a favourable \ntrade-off between representational richness, data efficiency, and interpretability, and substantially reduces \ncomputational and memory demands relative to fully trainable deep learning models  [42, 44] . Rather than \napproximating an optimal end-to -end decoder, HeteroRC provides a principled and tractable state- space \ntransformation aligned with known properties of neural dynamics, supporting time -resolved and hypothesis -\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n27 \n \ndriven analyses in low-data, high-variability settings typical of electrophysiological research. \nIn conclusion, the HeteroRC framework provides an interpretable and computationally efficient foundation \nfor time-resolved neural decoding. Its lightweight and interpretable architecture, together with its compatibility \nwith time -resolved decoding, makes it well -suited for a broad range of electrophysiological applications, \nincluding but not limited to studies of internally generated cognitive states and latent neural representations. To \nfacilitate reproducibility and adoption by the community, the full HeteroRC decoding and interpretation codes \nare made openly available as a documented, open -source software package on GitHub  \n(https://github.com/rl671/heterorc). We anticipate that this resource will enable systematic investigation of how \nneural information is dynamically represented over time and promote methodological advances in neural \ndecoding beyond phase-locked responses, with relevance for both neuroscience and brain– computer interface \nresearch.\n \nMethods \nSimulated datasets \nTo systematically evaluate decoding sensitivity to distinct classes of neural dynamics, we generated \nsynthetic datasets mimicking commonly studied electrophysiological  regimes, including stimulus -locked \nevoked responses, non-phase -locked oscillatory activity, ISPC , and aperiodic (1/ f) activity. Simulations were \ndesigned to isolate each regime while maintaining realistic signal-to-noise characteristics and spatial structure, \nallowing controlled comparisons between decoding models. All codes used in this study can be found on GitHub \n(https://github.com/rl671/heterorc). \nEach simulated dataset consisted of 30 independent “subjects”. For each subject, we generated two -class \nclassification data with 40 trials per class (80 trials total). Signals were simulated for a standard 32-channel EEG \nmontage, sampled at 100 Hz, over a time window from 0 to 8 00 ms. Channel labels and regional groupings \nfollowed the international 10–20 EEG system. Frontal channels were defined as electrodes with labels beginning \nwith “F,” whereas posterior channels were defined as electrodes with labels beginning with “P,” “O,” or “CP.” \nFor simulations of evoked responses, induced oscillatory activity, and aperiodic spectral modulations, task -\nrelated signals were injected exclusively into posterior channels, with frontal channels serving as controls. In \nthe ISPC condition, oscillatory signals were injected into both frontal and posterior channel groups, and class \ndifferences were implemented by modulating the strength of phase synchronisation between these regions. \nAll simulations incorporated realistic background noise composed of a mixture of spatially uncorrelated \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n28 \n \nwhite noise and temporally correlated pink (1/ f) noise. For each trial and channel, white Gaussian noise and \npink noise were combined with fixed weights (0.4 and 0.6, respectively) and scaled to yield unit-variance signals \nprior to task-related modulation. This background noise was present in all condi tions and was identical across \nclasses, ensuring that classification performance depended exclusively on the injected physiological signal of \ninterest. \nFive signal regimes were simulated independently.  Across all regimes, task -related modulations were \nconfined to a post -stimulus time window between approximately 200 and 600 ms, with modest trial -to-trial \ntemporal jitter to avoid perfectly aligned onsets and offsets.  Across all simulated regimes, class identity was \nencoded via a moderate difference in the target feature (e.g., amplitude, phase synchronisation, or aperiodic \nparameters) embedded in a physiological noise background comprising a mixture of 1/f pink and white noise. \nThe resulting macroscopic signal differences are shown in the univariate contrasts (Figure 2). \nTo simulate stimulus-locked evoked activity, class-discriminative signals were implemented as Gaussian-\nshaped amplitude deflections time-locked to stimulus onset. The evoked response peaked around 400 ms post-\nstimulus, with small trial-to -trial temporal (± ~25 ms) and amplitude jitter.\n Although this temporal jitter was \nintroduced to approximate physiological variance, the macroscopic amplitude  deflections remained largely \nconsistent across trials, ensuring that discriminative information was predominantly phase-locked to the event. \nInduced activity was simulated as transient oscillatory bursts with a Gaussian temporal envelope (duration \n≈ 400 ms) at a target frequency (default 10 Hz). Oscillatory phase was randomised independently on each trial, \nrendering the signal non-phase-locked at the sensor level. Class information was encoded solely in oscillatory \namplitude. Bursts had  modest trial -wise variability in frequency (± 0.5 Hz) and amplitude to approximate \nphysiological variability.\n To assess frequency generality, additional simulations were performed at lower and \nhigher carrier frequencies (5, 15, and 25 Hz), while all other parameters were held constant. \nTo model phase-based functional connectivity (i.e., ISPC), oscillatory bursts were simultaneously injected \ninto frontal and posterior channel groups. For one class, the phase difference between frontal and posterior \nsignals was tightly clustered (high phase synchroni sation), whereas for the other class phase differences were \nbroadly distributed. Importantly, marginal oscillatory power was matched across classes; discriminative \ninformation was carried exclusively by inter -regional phase consistency rather th an local amplitude or evoked \nresponses.\n As for induced activity, ISPC simulations were primarily conducted at 10 Hz, with additional carrier \nfrequencies (5, 15, and 25 Hz) examined in supplementary analyses to verify frequency-independent decoding \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n29 \n \nperformance. \nAperiodic neural activity was simulated using a physiologically motivated spectral rotation model. For \neach channel, broadband noise was filtered to produce 1/f -like spectra with controllable slope (exponent) and \nintercept (offset in log-power space). Class differences were implemented as changes in either the spectral slope \nor intercept, pivoting around a fixed frequency ( default 10 Hz). Aperiodic modulations were applied using a \nsmooth temporal mask. Outside this interval, all channels exhibited identical baseline aperiodic structure. \nMotor imagery dataset \nWe evaluated HeteroRC on the publicly available Graz motor imagery dataset from BCI Competition IV  \n(2008), data set 2a (https://www.bbci.de/competition/iv/#datasets) [51]. The dataset comprises EEG recordings \nfrom nine healthy participants performing a cue -based motor imagery task. Each participant completed two \nrecording sessions on separate days; one session was designated as training data with class labels provided, and \nthe other as evaluation data. \nParticipants performed four types of motor imagery: left-hand movement, right-hand movement, both feet \nmovement, and tongue movement. Each session consisted of six runs, and each run included 48 trials (12 trials \nper class), yielding 288 trials per session. At the beginning of each trial, a fixation cross appeared on the screen \naccompanied by a brief auditory warning tone. After 2 s of fixation, a visual cue in the form of an arrow pointing \neither to the left, right, down, or up indicated the motor imagery class and remained on the screen for 1.25 s. \nParticipants were instructed to continue the motor imagery task for 2.75 s until the fixation cross disappeared. \nNo feedback was provided during the task. \nEEG was recorded using 22 Ag/AgCl electrodes arranged according to the international 10 –20 system, \nwith inter-electrode distances of approximately 3.5 cm. Signals were recorded monopolarly with the left mastoid \nas reference and the right mastoid as ground. Data were sampled at 250 Hz and bandpass-filtered online between \n0.5 and 100 Hz, with an additional 50 Hz notch filter applied to suppress line noise.\n Three additional EOG \nchannels were recorded for artefact monitoring but were not used for decoding analyses. \nOffline preprocessing was intentionally kept minimal. EEG data were bandpass-filtered between 0.5 and \n30 Hz, downsampled to 100 Hz, and re -referenced to the common average. No additional art efact correction \nprocedures were applied. \nAttentional priority mapping dataset \nWe further evaluated HeteroRC on a publicly available EEG dataset investigating history-driven attentional \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n30 \n \npriority maps  (https://osf.io/v7yhc ), originally reported by Duncan, van Moorselaar [36] . In this study, 24 \nparticipants performed a visual search task designed to induce statistical learning of spatial target probabilities, \nleading to the implicit formation of an attentional priority map. One participant was excluded due to a technical \nfailure during data loading, resulting in a final sample of 23 participants included in the analyses. \nParticipants completed a variant of the additional singleton task in which targets appeared more frequently \nat one of four cardinal locations (up, down, left, or right) across blocks, while other locations were less probable. \nCritically, on half of the trials, a task -irrelevant but salient visual “ping” stimulus was presented during the \nintertrial interval, whereas no such stimulus was shown on the remaining trials. This design explicitly allows \ncomparison of ping and no-ping conditions while holding task demands constant. \nEEG was recorded from 64 scalp electrodes arranged according to the international 10– 10 system using a \nBioSemi ActiveTwo system.\n In the original study, EEG signals were re-referenced to the common average, \nhigh-pass filtered at 0.01 Hz.  Epochs containing pronounced EMG or muscle -related art efacts were then \nidentified and removed. Independent component analysis (ICA) was subsequently applied to remove ocular and \nother stereotypical artefacts, and trials containing eye movements or saccades were further excluded based on \nconcurrent eye-tracking data.  \nTo ensure consistency across datasets and decoding analyses, we applied additional preprocessing steps to \nthe provided data. Specifically, EEG signals were bandpass-filtered between 0.5 and 30 Hz and downsampled \nto 100 Hz. No further arte fact rejection or channel selection was performed beyond the preprocessing \nimplemented in the original dataset (see Duncan, van Moorselaar [36]  for details). No time-domain baseline \ncorrection was applied in the decoding analyses, in order to preserve neural activity preceding the ping stimulus. \nUnivariate analyses of simulated data \nTo characterise the signal properties of the simulated datasets independently of decoding, we performed \nstandard univariate analyses at the subject level and summari sed results across simulated subjects. For each \nsimulation mode, data were generated for multiple simulated subjects with identical generative parameters but \nindependent noise realisations. \nEvoked responses were computed by averaging epoched sensor -level signals across trials and across all \nelectrodes within each subject and condition, yielding a global mean time course. Power spectral density (PSD) \nwas estimated using Welch’s method over the entire time window and averaged across trials and electrodes \nwithin each subject, with power expressed in decibels.\n To quantify phase-based interactions in the simulated \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n31 \n \ndata, we computed  ISPC, also referred to as the phase-locking value (PLV), between predefined frontal and \nposterior electrode groups. ISPC was computed as the magnitude of the average unit phase -difference vector \nacross trials within a predefined temporal window (0.2–0.6 s): \nISPC =∣ 1\n𝑁𝑁� 𝑒𝑒𝑖𝑖Δ𝜙𝜙𝑛𝑛\n𝑁𝑁\n𝑛𝑛=1\n∣, \nwhere Δ𝜙𝜙𝑛𝑛 denotes the phase difference between frontal and posterior regions for trial 𝑛𝑛, evaluated at the \nmidpoint of the temporal window, and 𝑁𝑁 is the number of trials. ISPC values range from 0 (random phase \ndifferences across trials) to 1 (perfect phase alignment). In addition, per-trial phase differences between regions \nwere pooled across subjects to visualis e the distribution of phase offsets. For all univariate measures, subject-\nlevel results were averaged across simulated subjects, and variability was quantified using the standard error of \nthe mean.  \nTime-resolved decoding and cross-temporal generalisation using LDA and SVM \nAs baseline decoding approaches, we employed LDA and linear SVM, which are widely used in time-\nresolved EEG/MEG/LFPs decoding. To evaluate instantaneous signal decodability,  decoding was performed \nindependently at each discrete time point using a sliding-estimator approach, such that classifiers were trained \nand evaluated on sensor -level amplitude patterns across all sensors at each individual time sample. Prior to \nclassification, features were standardised using z-scoring based on the training data only.  For decoding based \non conventional time -frequency features (Figure S2), we extracted alpha-band power (8 –12 Hz) using two  \nclassic approaches. In the first, we computed spectral power via a Morlet wavelet transform (with cycles \ndynamically set to f/2), averaging the resulting estimates across 8–12 Hz. In the second, we applied an 8–12 Hz \nbandpass filter to the raw signals and extracted instantaneous power via a Hilbert transform (calculated as the \nsquared absolute value of the analytical signal). For both approaches, the resulting power time series were \nstandardised and submitted to the identical sliding-estimator LDA pipeline. \nFor simulated datasets and the attentional priority mapping dataset, time -resolved decoding accuracy was \nestimated using stratified 5-fold cross-validation. At each time point, classifiers were trained on the training \nfolds and evaluated on held-out data, yielding a decoding accuracy time course that was then averaged across \nfolds. \nFor the motor imagery dataset, decoding performance was evaluated using a fixed train– test split defined \nby recording sessions. Classifiers were trained on data from the training session and tested on data from a \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n32 \n \nseparate recording session collected on a different day, providing a stringent assessment of cross-session \ngeneralisation without cross-validation. \nTo characterise the temporal stability and transformation of neural representations, we performed cross-\ntemporal generalisation analyses for both linear decoders for the motor imagery dataset . In this approach, \nclassifiers were trained at a given time point and tested at all other time points, yielding a two -dimensional \nmatrix of decoding accuracy as a function of training time and testing time. Cross-temporal generalisation was \nevaluated using the same cross-session train-test split as time-resolved decoding, with classifiers trained on one \nsession and tested on the other. This analysis assesses the temporal generalisation  of neural representations \nacross time under strict cross-session generalisation constraints. \nHeterogeneous Reservoir Computing (HeteroRC) \nOverview. HeteroRC is a temporal decoding framework based on an echo  state network for extracting task -\nrelevant information from multichannel electrophysiological time series. The model projects neural signals into \na high-dimensional recurrent state space with heterogeneous intrinsic timescales, followed by a linear readout \ntrained to de code experimental variables. By combining recurrent dynamics with explicit temporal \nheterogeneity, HeteroRC is sensitive to both stimulus-locked, phase -consistent activity and non-phase-locked \nneural dynamics that unfold over longer or variable temporal scales. \nHeteroRC differs from conventional recurrent neural networks in three key respects: (i) both the input-to-\nreservoir weights and the recurrent weights within the reservoir are fixed and not optimised during training, (ii) \nnon-linear transformations arise solely from the intrinsic reservoir dynamics rather than from learned \nrepresentations, and  (iii) temporal heterogeneity is explicitly imposed by assigning distinct intrinsic time \nconstants to reservoir units . As a result, HeteroRC functions as a structured dimension expansion of the input \ntime series, enabling information carried by evoked responses as well as induced, non- phase-locked dynamics \nto be represented within a common state space and decoded using a linear readout. \nInput representation and normalisation. Let 𝐗𝐗 ∈ ℝ\n𝑁𝑁×𝐶𝐶×𝑇𝑇denote the epoched neural time-series data, where \n𝑁𝑁 is the number of trials, 𝐶𝐶 is the number of channels, and 𝑇𝑇 is the number of time points. For each trial 𝑛𝑛, \nthe multivariate time series 𝐱𝐱𝑛𝑛(𝑡𝑡) ∈ ℝ𝐶𝐶 is provided directly to the reservoir at each time point 𝑡𝑡. To ensure \nstrict separation between training and testing data, input normali sation (across all trials, channels, and time \npoints) was performed independently within each cross-validation fold. Specifically, input amplitudes were \nscaled using the 99th percentile of the absolute signal amplitude computed from the training data only, and the \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n33 \n \nsame global scaling factor was applied to the corresponding test data.  This normalisation  constrains input \nmagnitudes to a consistent dynamic range, ensuring that reservoir units operate within the sensitive (non-\nsaturated) regime of the tanh nonlinearity. \nReservoir architecture and dynamics. The reservoir consists of 𝑅𝑅  recurrent units with leaky -integrator \ndynamics. The state of the reservoir at time 𝑡𝑡  for trial 𝑛𝑛  is denoted by 𝐫𝐫𝑛𝑛(𝑡𝑡) ∈ ℝ𝑅𝑅 . Reservoir dynamics \nevolve according to the discrete-time update equation: \n𝐫𝐫𝑛𝑛(𝑡𝑡) = (1 − 𝜶𝜶) ⊙ 𝐫𝐫𝑛𝑛(𝑡𝑡 −1) + 𝜶𝜶 ⊙tanh(𝐖𝐖in𝐱𝐱𝑛𝑛(𝑡𝑡) + 𝐖𝐖res𝐫𝐫𝑛𝑛(𝑡𝑡 −1) + 𝐛𝐛), \nwhere 𝐖𝐖in ∈ ℝ𝑅𝑅×𝐶𝐶 is a fixed random input projection matrix, 𝐖𝐖res ∈ ℝ𝑅𝑅×𝑅𝑅 is a fixed sparse recurrent weight \nmatrix (10% connectivity density), and 𝐛𝐛 ∈ ℝ𝑅𝑅  is a bias vector. The operator  ⊙  denotes element-wise \nmultiplication, and 𝜶𝜶 ∈ ℝ𝑅𝑅 is a vector of unit -specific leak rates.  The input weights 𝐖𝐖in, recurrent weights \n𝐖𝐖res , and biases 𝐛𝐛  are randomly initialized and remain  fixed throughout training and testing. Recurrent \nweights are scaled to achieve a target spectral radius less than unity, ensuring echo-state stability. \nHeterogeneous intrinsic timescales. A defining feature of HeteroRC is the explicit introduction of \nheterogeneous intrinsic timescales across reservoir units. Each reservoir unit 𝑖𝑖 is assigned a time constant 𝜏𝜏𝑖𝑖, \ndrawn from a log-normal distribution and constrained to lie within a physiologically plausible range \n[𝜏𝜏min, 𝜏𝜏max] . Distribution parameters are chosen such that the mode of the distribution corresponds to a \nbiologically motivated timescale (e.g., 𝜏𝜏mode = 10 ms), while allowing a heavy tail that captures slower \ndynamics. The leak rate for each unit is derived from its time constant according to: \n𝛼𝛼𝑖𝑖= 1 − 𝑒𝑒𝑒𝑒𝑒𝑒�− 1\n𝑓𝑓𝑠𝑠𝜏𝜏𝑖𝑖\n� \nwhere 𝑓𝑓𝑠𝑠 is the sampling frequency. Under this formulation, each reservoir unit effectively acts as a low-pass \ntemporal filter with a cutoff frequency 𝑓𝑓𝑐𝑐= (2𝜋𝜋𝜏𝜏𝑖𝑖)−1 . The reservoir therefore implements a multiscale \ntemporal filter bank, enabling integration of input information over diverse temporal windows without explicit \nfrequency-domain decomposition. \nBidirectional temporal processing. To enhance sensitivity to temporally extended and temporally symmetric \npatterns, HeteroRC employs bidirectional temporal processing. Each trial is processed independently in both \nforward and backward temporal directions using identical reservoir parameters. Backward reservoir states are \nobtained by applying the reservoir to time-reversed input sequences. Forward and backward reservoir states are \ncombined at each time point using element-wise multiplication: \n𝐫𝐫\n𝑛𝑛bi(𝑡𝑡) = 𝐫𝐫𝑛𝑛→(𝑡𝑡) ⊙ 𝐫𝐫𝑛𝑛←(𝑡𝑡) \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n34 \n \nThis multiplicative fusion emphasi ses temporally consistent structure across directions while suppressing \ntransient, direction-specific noise, and preserves the fixed-weight nature of the reservoir. Alternative processing \nstrategies, including unidirectional reservoirs and bidirectional fusion via element -wise averaging, were also \nevaluated in supplementary analyses (Figure S3 ). Unless otherwise stated, all main analyses in the text use \nbidirectional processing with multiplicative fusion. \nLinear readout and time -resolved decoding. Decoding is performed using a linear classifier trained on the \nbidirectional reservoir states at each time point. Let 𝐳𝐳𝑛𝑛(𝑡𝑡) = 𝐫𝐫𝑛𝑛bi(𝑡𝑡) denote the final reservoir representation. \nClass scores are computed as:  𝐬𝐬𝑛𝑛(𝑡𝑡) = 𝐖𝐖out\n⊤ 𝐳𝐳𝑛𝑛(𝑡𝑡), where 𝐖𝐖out is estimated using ridge-regulari sed linear \nclassification. Time-resolved decoding performance was evaluated using stratified 5-fold cross -validation for \nsimulated datasets and the attentional priority mapping dataset, and using a fixed cross-session train –test split \nfor the motor imagery dataset. Readout weights were trained on the training data and evaluated on held-out data \nat each time point, yielding a decoding accuracy time course that was averaged across folds or across test \nsamples, respectively. T o improve robustness of time -resolved decoding, decision scores were temporally \nsmoothed using a Gaussian kernel (FWHM = 25 ms) prior to class prediction.  Cross-temporal generalisation \nanalyses were performed only for the motor imagery dataset and followed the same cross-session train–test split. \nReadout weights trained at a given time point in the training session were evaluated across all time points in the \nindependent test session, yielding training-by-testing time generalisation matrices.  In this decoding procedure, \nonly the readout weights 𝐖𝐖\nout are trained. All reservoir parameters were fixed and shared between training \nand testing data. Consequently, decoding performance reflects the expressive capacity of the reservoir dynamics \nrather than representational learning in the readout. \nHeteroRC parameter settings. In this study, unless otherwise stated, the following parameter settings were \nused for all HeteroRC analyses. For simulated datasets, the reservoir size was set to 𝑅𝑅= 350  units. For \nanalyses of real EEG datasets, including the motor imagery and attentional priority mapping datasets, the \nreservoir size was increased to 𝑅𝑅= 800 units to accommodate the higher dimensionality and variability of \nempirical data. Across all datasets, the recurrent weight matrix was scaled to a target spectral radius of 0.95. \nIntrinsic time constants were drawn from a log -normal distribution with mode 𝜏𝜏\nmode = 0.01 s and shape \nparameter 𝜎𝜎= 0.8, producing a heavy-tailed distribution of timescales. Time constants were constrained to lie \nwithin the range [𝜏𝜏min, 𝜏𝜏max] = [0. 002,0.08]s. Bidirectional temporal processing was enabled in all analyses, \nwith forward and backward reservoir states combined using element-wise multiplication. As a general guideline \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n35 \n \nfor hyperparameter selection in this framework, 𝑅𝑅 should be chosen to be large enough to provide a sufficiently \nrich high-dimensional space for disentangling complex dynamics.  A common heuristic is to set 𝑅𝑅  to \napproximately 10 to 15 times the number of input channels. Additionally, setting the spectral radius just below \nunity (typically between 0.90 and 0.99) is a typical  practice to rigorously enforce the echo state property, \nensuring the network maintains a stable, fading memory of past inputs without transitioning into chaotic \ndynamics. \nIndividual-level interpretability framework \nTo interpret task -relevant information represented within the HeteroRC reservoir, we first developed a n \nindividual-level interpretability framework that progressively links reservoir dynamics to sensor -level \nneurophysiological signals. This framework extends covariance -based activation pattern analysis to recurrent \nreservoirs and enables separation of heterogeneous temporal,  spectral, and spatial contributions to decoding \nperformance. \nIdentification of task-relevant latent dynamics in reservoir space. We first identified reservoir units that \ncontributed most strongly to the decoding decision. In multivariate decoding models, linear readout weights \ncannot be directly interpreted as feature importance because they reflect both signal encoding and noise \nsuppression induced by feature covariance [18]. To recover task-related activation patterns within the reservoir, \nwe applied a covariance-based transformation of the classifier weights, following the approach introduced by \nHaufe and colleagues [18] to yield interpretable activation coefficients. \nLet 𝐫𝐫(𝑡𝑡) ∈ ℝ\n𝑅𝑅  denote the reservoir state vector at time 𝑡𝑡 , and let 𝑦𝑦 �(𝑡𝑡) = 𝐖𝐖out\n⊤ 𝐫𝐫(𝑡𝑡)  denote the \ncorresponding decision value of the linear readout. Activation patterns 𝐀𝐀 were computed as: \n𝐀𝐀= 𝚺𝚺𝐫𝐫𝐖𝐖out𝚺𝚺𝑦𝑦�\n−1 \nwhere 𝚺𝚺𝐫𝐫 is the covariance matrix of reservoir states and 𝚺𝚺𝑦𝑦�  denotes the covariance of the decision values. \nThis transformation yields activation coefficients reflecting how strongly each reservoir unit encodes task -\nrelevant variance. \nFor multi-class decoding, a global importance score was assigned to each reservoir unit as the maximum \nabsolute activation coefficient across classes. The dominant class contribution for each unit was defined as the \nclass associated with the largest absolute activation coefficient. To disentangle distinct dynamical mechanisms \nrepresented within the high-dimensional reservoir state space, we selected the most informative units and \nperformed functional clustering based on their temporal activity profiles.  For individual-level interpretations, \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n36 \n \nwe retained the top 25 units by importance. Prior to clustering, the time series of the selected units were z-score \nstandardised to ensure that grouping was driven by the shape of temporal dynamics rather than signal amplitude. \nWe then applied hierarchical clustering using Ward’s method with the Euclidean distance metric. This procedure \ngrouped reservoir units into a small number of functional clusters reflecting distinct temporal and spectral \ndynamics. For each cluster 𝑘𝑘, a virtual source signal 𝑣𝑣𝑘𝑘(𝑡𝑡) was constructed as a weighted combination of each \nconstituent reservoir unit 𝑖𝑖 within cluster 𝒞𝒞𝑘𝑘: \n𝑣𝑣𝑘𝑘(𝑡𝑡) = 1\n∣ 𝒞𝒞𝑘𝑘∣ � 𝑤𝑤𝑖𝑖\n𝑖𝑖∈𝒞𝒞𝑘𝑘\n sign(𝑎𝑎𝑖𝑖,dom) 𝑟𝑟𝑖𝑖(𝑡𝑡), \nwhere 𝑤𝑤𝑖𝑖 denotes the unit importance score and the sign term enforces consistency with the dominant class \ncontribution. These virtual sources represent low-dimensional summaries of task-relevant reservoir dynamics. \nFor each virtual source 𝑣𝑣𝑘𝑘(𝑡𝑡), we quantified its task -relevant dynamics in the time, time –frequency, and \nspectral domains. First, event-related responses were computed by averaging 𝑣𝑣𝑘𝑘(𝑡𝑡) within each experimental \ncondition across trials, yielding class-wise virtual -source evoked responses. Second, time –frequency \nrepresentations were estimated using complex Morlet wavelets over frequencies from 2 to 40 Hz (2  Hz steps). \nThe number of cycles scaled with frequency (𝑛𝑛cycles = 𝑓𝑓/2), and power was computed for each condition and \nthen averaged across trials. To facilitate interpretation of induced changes, power was baseline-corrected in \ndecibel units by referencing each frequency bin to the mean power in a pre-stimulus baseline window (i.e., \n10log 10(𝑃𝑃/𝑃𝑃base) ). Third, PSD was computed from 𝑣𝑣𝑘𝑘(𝑡𝑡)  using Welch’s method and condition-specific \nPSDs were averaged across trials and visuali sed on a log-frequency axis. To dissociate aperiodic and periodic \ncomponents, spectra were parameterised using FOOOF [28] over 2–40 Hz with a maximum of three oscillatory \npeaks, and the fixed aperiodic mode was used to summarise 1/f-like structure. \nBack-projection of reservoir dynamics to sensor space. In this stage, we linked task -relevant reservoir \ndynamics back to the sensor space to recover their spatial origins . Under the linear forward model of volume \nconduction, the spatial topography of a latent source can be estimated by its covariance with sensor -level \nmeasurements [18, 59, 60]. To capture the spatial expression of the reservoir dynamics specifically during the \ndecision-making process, we restricted our analysis to a focused time window ( ±100 ms) centred around the \npeak decoding accuracy . We estimated a spatial projection 𝐩𝐩𝑘𝑘∈ ℝ𝐶𝐶  for each virtual source 𝑣𝑣𝑘𝑘(𝑡𝑡)  by \ncomputing its covariance with the raw sensor signals 𝐱𝐱(𝑡𝑡) within this window: \n𝐩𝐩𝑘𝑘= Cov(𝐱𝐱(𝑡𝑡), 𝑣𝑣𝑘𝑘(𝑡𝑡)), \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n37 \n \nThe resulting spatial patterns 𝐩𝐩𝑘𝑘 represent the topographical distribution of the specific latent dynamics (i.e., \nthe virtual source signal) corresponding to each functional cluster.  \nGroup-level interpretability framework \nTo determine whether the latent representational dynamics observed at the single -subject level were \nconserved across the cohort, we extended the interpretability framework to perform group -level inference. \nBecause randomly initialised reservoirs yield specific and incommensurable latent spaces, direct cross-subject \nalignment of reservoir units is challenging. To overcome this, we developed a sensor -space spatial matching \napproach, leveraging the fact that the physical EEG sensor array is anatomically standardised across participants. \nFor each participant, we extracted the 10 most informative reservoir units and computed their spatial \ntopographies (sensor-level covariance maps) within a ±100ms window centred on the participant-specific peak \ndecoding time. As two reservoir units may capture identical neural dynamics but with inverted polarities, we \napplied a spatial sign -alignment procedure: for each unit’s topography, we identified the sensor with the \nmaximum absolute covariance and multiplied the entire topographical map by -1 if this peak value was negative. \nThe corresponding sign-flip coefficient was retained for subsequent signal reconstruction. These sign -aligned \nmaps were then z-scored to ensure that matching was driven strictly by the underlying spatial distribution (i.e., \ncortical generator patterns) rather than absolute covariance amplitude.  Standardised spatial maps from all \nselected units across all participants were concatenated into a single global feature matrix. We then applied \nhierarchical clustering (Ward’s method, Euclidean distance) to group the units into distinct global spatial clusters. \nTo reconstruct the group-level neurophysiological dynamics, we mapped the globally clustered units back \nto their respective participants. Participant-specific virtual sources were reconstructed by computing a weighted \nsum of the constituent unit time series. During this reconstruction, each unit's time series was multiplied by its \npreviously determined sign-flip coefficient. Subject-level temporal, spectral, and time –frequency \nrepresentations were then extracted from these aligned virtual sources and averaged across participants to yield \ngrand-average profiles. Finally, to generate the grand-average spatial topography for each global cluster, the raw \nspatial covariance maps of the aligned virtual sources were averaged across participants and z-scored. \nTime-resolved decoding using RNN, LSTM, transformer, and EEGNet \nTo compare HeteroRC against fully trainable ANNs, we implemented time-resolved decoding pipelines \nbased on RNNs, LSTM networks, Transformers, and EEGNet. We restricted this comparison  to the simulated \nphase-locked evoked responses  because these highly consistent signals represent the most straightforward \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n38 \n \nscenario for ANNs to achieve stable performance. The data generation and cross-validation procedures \n(stratified 5-fold) were identical to those described above.  To preserve strict train/test separation, input \namplitudes were normalised independently within each fold using the 99th percentile of the absolute amplitude \ncomputed exclusively on the training data, matching the HeteroRC preprocessing pipeline. \nFor the RNN, LSTM, and Transformer, we evaluated two decoding formulations. In the standard \n(sequence-to-sequence) formulation, the full trial sequence was provided as input, and the models were trained \nto produce a continuous trajectory of class predictions, generating a discrete decision at every individual time \nstep based on the accumulated signal history. In the windowed (sequence-to -one) formulation, decoding was \nperformed using sliding windows of 10 samples (100 ms). To enable evaluation at early time points, inputs were \nzero-padded at the epoch onset. At each time step, windowed models output a single prediction based on the \npreceding 100 ms context. EEGNet was evaluated exclusively using this windowed formulation to respect its \nconvolutional architecture. \nRegarding model architectures, the standard RNN and LSTM models consisted of a single recurrent layer \n(with tanh nonlinearity for the RNN) followed by dropout and a linear classification head. The standard \nTransformer projected sensor-level inputs to a latent embedding, applied sinusoidal positional encoding, and \nprocessed the sequence through a two -layer Transformer encoder prior to linear classification. The windowed \nvariants retained these core architectures but aggregated temporal information to produce a single prediction \nper window: the windowed RNN and LSTM utili sed the final hidden state, while the windowed T ransformer \napplied mean pooling across the encoded time dimension. Hidden dimensionality was fixed at 32 units across \nall RNN, LSTM, and Transformer models to match representational capacity. \nAll models were implemented in PyTorch\n and trained on a NVIDIA GeForce RTX 4090 GPU. The RNN, \nLSTM, and Transformer models were trained using the Adam optimiser (learning rate = 0.002, weight decay = \n0.0001) for 30 training iterations with a dropout rate of 0.5. The batch size was set to 64 for standard models \nand 128 for windowed models. EEGNet was trained using the AdamW optimiser (learning rate = 0.005, weight \ndecay = 0.001) for 25 training iterations with a batch size of 64. Time-resolved decoding accuracy was computed \nindependently at each time point by comparing predicted and ground-truth class labels, yielding a time-resolved \naccuracy profile mirroring the evaluation approach used for HeteroRC. \nStatistical analysis \nStatistical significance of time-resolved decoding accuracy and cross-temporal generalisation matrices was \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n39 \n \nassessed using nonparametric cluster -based permutation tests across subjects, as implemented for one-sample \ntests against chance level. This approach controls for multiple comparisons across time points while making \nminimal assumptions about the underlying distributions. \nComparisons of peak decoding accuracy across decoding methods were performed using paired-sample t-\ntests (two-tailed), applied separately to predefined temporal windows and conditions. Unless otherwise stated, \nstatistical tests were conducted at a significance threshold of p < 0.05. \nAcknowledgements \nThis project was supported by UKRI MRC intramural funding MC_UU_00030/15 to A.W. R.L. was supported \nby a Gates Cambridge Scholarship (OPP1144) and a postdoctoral fellowship from the Canadian Institutes for \nHealth Research (200883). S.L. was supported by Vetenskapsrådet under award 2023-00493. For the purpose \nof open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted \nManuscript version arising from this submission. \nAuthor contributions \nR.L.: Conceptualisation, Methodology, Software , Formal analysis, Data Curation , Writing - Original Draft , \nWriting - Review & Editing, Visualisation, Project administration. S.L.: Software, Validation, Formal analysis, \nWriting - Review & Editing. Y. L.: Software, Writing - Review & Editing. J.D.: Writing - Review & Editing. \nR.N.H.: Methodology, Writing - Review & Editing. A.W.: Methodology, Writing - Review & Editing, \nSupervision, Funding acquisition. \nDeclaration of interests \nThe authors declare no competing interests. \nReferences \n1. Haxby, J.V ., A.C. Connolly, and J.S. Guntupalli, Decoding neural representational spaces using multivariate \npattern analysis. Annu Rev Neurosci, 2014. 37: p. 435-56. \n2. Norman, K.A., et al., Beyond mind- reading: multi-voxel pattern analysis of fMRI data. Trends Cogn Sci, 2006. \n10(9): p. 424-30. \n3. Haxby, J.V ., et al., Distributed and overlapping representations of faces and objects in ventral temporal cortex. \nScience, 2001. 293(5539): p. 2425-30. \n4. Kamitani, Y . and F. Tong, Decoding the visual and subjective contents of the human brain.  Nat Neurosci, 2005. \n8(5): p. 679-85. \n5. Carlson, T.A., P. Schrater, and S. He, Patterns of activity in the categorical representations of objects.  J Cogn \nNeurosci, 2003. 15(5): p. 704-17. \n6. Grootswagers, T., S.G. Wardle, and T.A. Carlson, Decoding Dynamic Brain Patterns from Evoked Responses: A \nTutorial on Multivariate Pattern Analysis Applied to Time Series Neuroimaging Data. J Cogn Neurosci, 2017. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n40 \n \n29(4): p. 677-697. \n7. Lotte, F., et al., A review of classification algorithms for EEG-based brain-computer interfaces: a 10 year update. \nJ Neural Eng, 2018. 15(3): p. 031005. \n8. Ding, Y ., et al., EEG-based brain-computer interface enables real-time robotic hand control at individual finger \nlevel. Nat Commun, 2025. 16(1): p. 5401. \n9. King, J.R. and S. Dehaene, Characterizing the dynamics of mental representations: the temporal generalization \nmethod. Trends Cogn Sci, 2014. 18(4): p. 203-10. \n10. Peelen, M.V . and P.E. Downing, Testing cognitive theories with multivariate pattern analysis of neuroimaging \ndata. Nat Hum Behav, 2023. 7(9): p. 1430-1441. \n11. Robinson, A.K., G.L. Quek, and T.A. Carlson, Visual Representations: Insights from Neural Decoding. Annu Rev \nVis Sci, 2023. \n12. Cichy, R.M., D. Pantazis, and A. Oliva, Resolving human object recognition in space and time.  Nat Neurosci, \n2014. 17(3): p. 455-62. \n13. Hebart, M.N. and C.I. Baker, Deconstructing multivariate decoding for the study of brain function. Neuroimage, \n2018. 180(Pt A): p. 4-18. \n14. Bae, G.Y . and S.J. Luck, Dissociable Decoding of Spatial Attention and Working Memory from EEG Oscillations \nand Sustained Potentials. J Neurosci, 2018. 38(2): p. 409-422. \n15. Lu, R., et al., Parietal alpha stimulation causally enhances attentional information coding in evoked and \noscillatory activity. Brain Stimul, 2025. 18: p. 114-27. \n16. Renton, A.I., D.R. Painter, and J.B. Mattingley, Optimising the classification of feature -based attention in \nfrequency-tagged electroencephalography data. Scientific Data, 2022. 9(1). \n17. Trammel, T., et al., Decoding semantic relatedness and prediction from EEG: A classification method comparison. \nNeuroimage, 2023. 277: p. 120268. \n18. Haufe, S., et al., On the interpretation of weight vectors of linear models in multivariate neuroimaging.  \nNeuroimage, 2014. 87: p. 96-110. \n19. Pantazis, D., et al., Decoding the orientation of contrast edges from MEG evoked and induced responses.  \nNeuroimage, 2018. 180(Pt A): p. 267-279. \n20. Foster, J.J. and E. Awh, The role of alpha oscillations in spatial attention: limited evidence for a suppression \naccount. Curr Opin Psychol, 2019. 29: p. 34-40. \n21. Stecher, R., R.M. Cichy, and D. Kaiser, Decoding the rhythmic representation and communication of visual \ncontents. Trends Neurosci, 2025. 48(3): p. 178-188. \n22. Lundqvist, M., et al., Beta: bursts of cognition. Trends Cogn Sci, 2024. 28(7): p. 662-676. \n23. Johnson, E.L., et al., A rapid theta network mechanism for flexible information encoding.  Nat Commun, 2023. \n14(1): p. 2872. \n24. Siegel, M., T.H. Donner, and A.K. Engel, Spectral fingerprints of large -scale neuronal interactions.  Nat Rev \nNeurosci, 2012. 13(2): p. 121-34. \n25. Fries, P., Rhythms for Cognition: Communication through Coherence. Neuron, 2015. 88(1): p. 220-35. \n26. Palva, J.M., et al., Neuronal synchrony reveals working memory networks and predicts individual memory \ncapacity. Proc Natl Acad Sci U S A, 2010. 107(16): p. 7580-5. \n27. Albouy, P., et al., Supramodality of neural entrainment: Rhythmic visual stimulation causally enhances auditory \nworking memory performance. Sci Adv, 2022. 8(8): p. eabj9782. \n28. Donoghue, T., et al., Parameterizing neural power spectra into periodic and aperiodic components. Nat Neurosci, \n2020. 23(12): p. 1655-1665. \n29. Lu, R., et al., Aperiodic and oscillatory systems underpinning human domain-general cognition. Communications \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n41 \n \nBiology, 2024. 7: p. 1643. \n30. Lu, R., E. Pollitt, and A. Woolgar, Distinct and complementary mechanisms of oscillatory and aperiodic alpha \nactivity in visuospatial attention. Imaging Neuroscience, 2025. 3. \n31. Wolff, M.J., et al., Dynamic hidden states underlying working -memory-guided behavior. Nat Neurosci, 2017. \n20(6): p. 864-871. \n32. Barbosa, J., et al., Interplay between persistent activity and activity -silent dynamics in the prefrontal cortex \nunderlies serial biases in working memory. Nat Neurosci, 2020. 23(8): p. 1016-1024. \n33. Trubutschek, D., et al., Probing the limits of activity-silent non-conscious working memory. Proc Natl Acad Sci U \nS A, 2019. 116(28): p. 14358-14367. \n34. Stokes, M.G., 'Activity-silent' working memory in prefrontal cortex: a dynamic coding framework.  Trends Cogn \nSci, 2015. 19(7): p. 394-405. \n35. Rose, N.S., et al., Reactivation of latent working memories with transcranial magnetic stimulation. Science, 2016. \n354(6316): p. 1136-39. \n36. Duncan, D.H., D. van Moorselaar, and J. Theeuwes, Pinging the brain to reveal the hidden attentional priority \nmap using encephalography. Nat Commun, 2023. 14(1): p. 4749. \n37. Barbosa, J., D. Lozano -Soldevilla, and A. Compte, Pinging the brain with visual impulses reveals electrically \nactive, not activity-silent, working memories. PLoS Biol, 2021. 19(10): p. e3001436. \n38. Karimi-Rouzbahani, H., et al., Temporal Variabilities Provide Additional Category-Related Information in Object \nCategory Decoding: A Systematic Comparison of Informative EEG Features.  Neural Comput, 2021. 33 (11): p. \n3027-3072. \n39. Karimi-Rouzbahani, H. and A. Woolgar, When the Whole Is Less Than the Sum of Its Parts: Maximum Object \nCategory Information and Behavioral Prediction in Multiscale Activation Patterns. Front Neurosci, 2022. 16: p. \n825746. \n40. Roy, Y ., et al., Deep learning-based electroencephalography analysis: a systematic review. J Neural Eng, 2019. \n16(5): p. 051001. \n41. Varoquaux, G., Cross-validation failure: Small sample sizes lead to large error bars. Neuroimage, 2018. 180(Pt \nA): p. 68-77. \n42. Jaeger, H., The “echo state” approach to analysing and training recurrent neural networks, in German national \nresearch center for information technology gmd technical report. 2001: Bonn, Germany. \n43. Maass, W., T. Natschläger, and H. Markram, Real-Time Computing Without Stable States: A New Framework for \nNeural Computation Based on Perturbations. Neural Comput, 2002. \n44. Yan, M., et al., Emerging opportunities and challenges for the future of reservoir computing. Nat Commun, 2024. \n15(1): p. 2056. \n45. Verstraeten, D., et al. The unified Reservoir Computing concept and its digital hardware implementations. in 2006 \nEPFL LATSIS Symposium. 2006. EPFL, Lausanne. \n46. Suarez, L.E., et al., Connectome-based reservoir computing with the conn2res toolbox. Nat Commun, 2024. 15(1): \np. 656. \n47. Li, G., S. Li, and X.J. Wang, A hierarchy of time constants and reliable signal propagation in the marmoset \ncerebral cortex. Nat Commun, 2025. 16(1): p. 11640. \n48. Dahmen, D., et al., How heterogeneity shapes dynamics and computation in the brain. Neuron, 2025. \n49. Spitmaan, M., et al., Multiple timescales of neural dynamics and integration of task-relevant signals across cortex. \nProc Natl Acad Sci U S A, 2020. 117(36): p. 22522-22531. \n50. Bernacchia, A., et al., A reservoir of time constants for memory traces in cortical neurons.  Nat Neurosci, 2011. \n14(3): p. 366-72. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n42 \n \n51. Tangermann, M., et al., Review of the BCI Competition IV . Front Neurosci, 2012. 6: p. 55. \n52. Grossberg, S., Recurrent neural networks. Scholarpedia, 2013. 8: p. 1888. \n53. Graves, A., Long Short-Term Memory, in Supervised Sequence Labelling with Recurrent Neural Networks , A. \nGraves, Editor. 2012, Springer Berlin Heidelberg: Berlin, Heidelberg. p. 37-45. \n54. Ashish Vaswani, et al., Attention Is All You Need, in Advances in Neural Information Processing Systems. 2017. \n55. Lawhern, V .J., et al., EEGNet: a compact convolutional neural network for EEG-based brain–computer interfaces. \nJournal of Neural Engineering, 2018. 15(5): p. 056013. \n56. Cohen, M.X., Analyzing Neural Time Series Data: Theory and Practice. 2014, The MIT Press. \n57. Jeon, Y ., et al., Event -related (De)synchronization (ERD/ERS) during motor imagery tasks: Implications for \nbrain–computer interfaces. International Journal of Industrial Ergonomics, 2011. 41(5): p. 428-436. \n58. Song, Y ., et al., EEG Conformer: Convolutional Transformer for EEG Decoding and Visualization. IEEE Trans \nNeural Syst Rehabil Eng, 2023. 31: p. 710-719. \n59. Hämäläinen, M., et al., Magnetoencephalography —theory, instrumentation, and applications to noninvasive \nstudies of the working human brain. Reviews of Modern Physics, 1993. 65(2): p. 413-497. \n60. Hauk, O., M. Stenroos, and M.S. Treder, Towards an objective evaluation of EEG/MEG source estimation methods \n- The linear approach. Neuroimage, 2022. 255: p. 119177. \n \n  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n43 \n \nSupplementary \n \n \nFigure S1. HeteroRC decodes non-phase-locked neural dynamics across a broad range of frequencies. \nTime-resolved decoding accuracy for simulated induced oscillatory power (left) and inter -site phase clustering \n(ISPC; right) when task -relevant modulations were ce ntred at frequencies (5 Hz, 15 Hz, and 25 Hz) different \nfrom those used in the main text.  HeteroRC robustly decodes task -relevant information across frequencies, \nduring the time  window that they were introduced (0.2-0.6s), whereas linear raw amplitude-based decoders \nremain at chance. This demonstrates that HeteroRC decoding performance does not rely on tuning to a specific \noscillatory band. Conventions as Figure 2. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n44 \n \n \nFigure S2. Conventional time -frequency feature extraction restricts decoding sensitivity and induces \ntemporal smearing.  \nComparison of linear discriminant analysis (LDA) performance applied to a priori narrow-band power extracted via \nMorlet wavelets (8–12 Hz; orange) and a Hilbert filter approach (8–12 Hz; green). Panels display decoding accuracy \nfor simulated data in which the two classes varied in their phase-locked evoked responses, induced oscillatory power \ncentred at 10 Hz and 20 Hz, inter-site phase clustering (ISPC) at 10 Hz, and aperiodic spectral modulations affecting \nthe 1/f slope and intercept  (as in Figure 2) . Simulated task-relevant modulations were confined to a 0.2 –0.6 s time \nwindow. Conventions as in Figure 2. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n45 \n \n \nFigure S3. Bidirectional reservoir processing improves temporal precision and reduces decoding smear. \nComparison of unidirectional processing  (left column)  and bidirectional processing  (right column)  with state \naveraging for decoding of phase-locked evoked responses (a), induced oscillatory power with randomized phase (b), \ninter-site phase clustering (ISPC) (c), and aperiodic spectral modulations affecting the 1/f slope (d) and offset (e) . \nSimulated task-relevant modulations were confined to a 0.2 –0.6 s time window. Unidirectional processing exhibits \nsystematic temporal lag and smearing in decoding peaks, whereas bidirectional processing improves temporal \nalignment. Yet, averaging-based fusion shows residual temporal smearing for aperiodic modulations  (d and e) . In \ncontrast, bidirectional multiplicative fusion (used in the main text, Figure 2 ) does not exhibit this smearing.  \nConventions as in Figure 2.  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n46 \n \n \n \n \nFigure S4. Individual-subject decoding results for the motor imagery dataset.  \nTime-resolved decoding accuracy for each participant in the motor imagery dataset, shown for HeteroRC  (red) and \nlinear discriminant analysis (LDA, blue). For visualisation purposes, decoding curves were smoothed with a Gaussian \nkernel (50 ms FWHM). These plots illustrate the consistency and inter -subject variability of decoding dynamics \nacross participants.  \n \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint \n\n \n47 \n \n \nFigure S5. Individual-subject decoding results for the attentional priority mapping task. \nTime-resolved decoding accuracy for each participant in the attentional priority mapping dataset, shown separately \nfor Ping and No -Ping conditions and for HeteroRC and LDA. Decoding curves were smoothed with a Gaussian \nkernel (50 ms FWHM) for visualisation.  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted April 7, 2026. ; https://doi.org/10.64898/2026.04.04.716475doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}